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PDES 


Set Books 


G. D. Smith, Numerical Solution of Partial Differential Equations (Oxford, 1971). 
H. F. Weinberger. A First Course in Partial Differential Equations (Blaisdell, 1965). 


It is essential to have these books: the course is based on them and will not make sense 
without them. They are referred to in the text as S and W respectively. 

Unit 5 is based on S: Chapter 1, pages 6 to 8 and Chapter 2. pages 10 to 24, 32 to 
40 and 46 to 52. 


Conventions 


Before working through this text make sure you have read A Guide to the Course: 
Partial Differential Equations of Applied Mathematics. References to Open University 
courses in mathematics take the form: 


Unit M100 13, Integration II for the Mathematics Foundation Course, 
Unit M201 23, The Wave Equation for the Linear Mathematics Course. 
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5.0 INTRODUCTION 


The various analytical methods of obtaining solutions to partial differential equations, 
some of which are discussed in this course, can be applied only to a restricted range of 
problems. For example, the initial value problem 


ди ?u . 
gr 9 ga À + (sin xu, t) = 0 0<x<l,t>0 


u(0, t) = u(1,2) = 0 t>0 
u(x, 0) = f(x) 0x 


cannot be solved by the method of separation of variables or by means of integral 
transforms (e.g. Laplace and Fourier transforms). Even when solutions can be found 
analytically, difficulties often arise in evaluating them numerically. For instance, 
the solution may involve finding the Fourier coefficients of the function which appears 
in an initial condition (Unit M201 32, The Heat Conduction Equation, Section 32.2). 
Each coefficient involves an integral which implies that, if we are aiming for full 
numerical precision, an infinitude of integrals must be evaluated. Unless the given 
function is particularly simple (for example, a polynomial, exponential, or trigono- 
metric function) only a finite number of these integrations can be carried out in a 
finite time. We are therefore reduced to approximating the Fourier series (e.g. by 
truncating it), and might even have to estimate approximately each of the separate 
integrals in the series. As a result we would end with only an approximate solution 
to our original problem. 


x & l, 


In this course we present some numerical methods which can be applied to a wide class 
of problems, but which inherently give only approximate values for the solutions. 
We do not claim that the numerical methods we shall illustrate can be applied to all 
problems involving partial differential equations or that they are the only methods 
available, but they are currently the methods most frequently used and most widely 
applicable. 


In Unit M201 7, Recurrence Relations and Unit M20! 8, Numerical Solution of 
Simultaneous Algebraic Equations we saw that, in numerical mathematics, attention 
had to be paid on the one hand to the formulation of the problem, and on the other 
hand to the method of solution. 


We have to consider the possibility of inherent instability in the problem, in which 
case small changes in the data make large changes in the solution. We have met this 
previously in this course, and we described such a problem as not being continuous 
with respect to its data. 


Throughout this unit we shall concentrate on initial value problems, i.e. problems in 
which the value of the solution (or its normal derivative) is given along the initial line 
(t = 0). Initial value problems in which a boundary condition is specified for t > 0 
at each point of the spatial boundary are properly posed for parabolic equations. 
In particular they exhibit continuity, and so do not suffer from inherent instability, 
as we have seen in Unit 3, Elliptic and Parabolic Equations. Similarly, the Cauchy 
problem (i.e. the problem as above, but with the solution and its normal derivative 
both specified on the initial line) is properly posed for hyperbolic equations. (We 
usually refer to these problems as initial-bowidary value problems.) By analogy with 
the treatment of initial-value problems for ordinary differential equations, discussed 
in Unit M201 21, Numerical Solution of Differential Equations we use approximating 
recurrence relations ina step-by-step process leading to a simple system of simultaneous 
algebraic equations. Each of the unknowns in the system of equations represents an 
approximate value of the solution to the differential equation at some predetermined 
point. We might expect that as the number of points is increased the numerical solution 
becomes more accurate. We shall employ methods for which we can prove (under 
suitable conditions) that this is the case. 


The recurrence relations are obtained by finite-difference methods, and we shall discuss 
explicit and implicit schemes. In the former, as the name implies, each equation 
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expresses one unknown value in terms of known values. In the latter, however, several 
unknown values must be solved for simultaneeusly. 


Following the discussion of the numerical solution of ordinary differential equations 
in Unit M201 24, we might expect that some of our methods exhibit induced instability. 
In this unit we shall illustrate such instability with specific numerical examples, and 
as a result we shall gain some insight into the factors which induce this phenomenon. 
We shall support this insight by a full theoretical treatment in Unit 8, Stability. 


In the main part of this unit we shall produce and use only the simplest finite-difference 
equations for parabolic and hyperbolic equations and, if necessary, for the initial and 
boundary conditions. More accurate schemes can be constructed, but we shall not 
discuss them here. Their existence illustrates how the numerical analyst is constantly 
searching for methods of greater efficiency and wider applicability. 


Numerical analysis has grown rapidly over the past few decades owing to the increasing 
availability of high-speed digital computers. The methods we shall investigate are 
designed primarily for such machines and in some SAQs in this unit you will be 
asked to run simple computer programs. These programs have already been written 
and are available as library routines; and the investigation of the results of computer 
runs using these programs is an integral part of the course. Occasionally. you may also 
be asked to write small programs in BASIC, but we shall not make it mandatory 
for you to do so. 
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5.1 ELEMENTARY FINITE-DIFFERENCE APPROXIMATIONS 


We begin our work in numerical analysis by revising some ideas which form the 
basis of the finite-difference method as applied to partial differential equations. We 
shall need to use Taylor series (Unit M 100 14, Sequences and Limits I and Unit M201 
14, Bilinear and Quadratic Forms) and finite differences (Unit M100 4, Finite Dif- 
ferences). 


To see how these ideas are used in practice we give a broad outline of the finite- 
difference method before going into details. The method requires us, at first, to 
choose a finite set of values for each of the domain variables in the problem. We 
usually choose these values in a systematic way in order to simplify the resulting 
theory and practice. For example, if we consider a problem with just two independent 
variables, denoted by x and t, then we take the sets of values (x;:i —0, L, .. . and 
ít 1) = 0,1,...}, where successive members of each set are related by the formulas 


Хы =X, +h 
and 

баа = 0+ К 
Here h апа К are constants which we can choose at our discretion, subject to certain 
conditions, as we shall see. хо and го are chosen so that the set of points 


(0:20) = (xo + ih, to + JK} 


"covers" a suitable domain. In practice, we frequently choose coordinates such that 
(Xo. to) = (0,0). We shall sometimes denote the point (x;, tj) simply as i,j where 
Хо, to, h and k are understood. We refer to the set of points {(x;, t,)} as a mesh or 
grid, and we call h and k the mesh lengths in the x- and t-directions respectively. 
The values u(x;, 1) of a function и at the mesh points are called pivotal values. 


Since it is possible to perform only a finite number of calculations in a finite time (even 
by computer!) we have to restrict the solution domain for any initial value problem 
that we care to tackle. From now on we shall assume that the mesh is finite; say 
{(x;. t:i = 0,12,..., N and j = 0,1.2,..., М). 


It is always possible to measure x in suitable units so that x, = 1. (This is discussed 
in S: pages 9 and 10, which is not, however, set as a reading passage.) Then 


h2— 
N 
where N + 1 is the number of mesh points in the x-direction. The solution domain 
for the finite-difference calculations for two independent variables x and t in a 
rectangular region is shown. 


Mk ae 


k 
3k h 
2k 
k 
> 
0 h 2h. 3M mE -Nh x 
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The partial derivatives in the partial differential equation to besolved are approximated 
at each point by finite-dillerence formulas based on the chosen mesh and obtained 
using Taylor approximations. You have seen these ideas before, applied to ordinary 
differential equations in Section 21.2 of Unit M20/ 21. 


Taylor's Theorem (Unit M /00 14) is that 
Г 
u(x + Л) = u(x) + hu(x) +... + = ми) + CU) 
where 


IC.) < c4 B, ibt "|, 
provided 

Ш" 909) < Ba, Felxx +h] 
and «"*" is continuous throughout [x, x + h]—or [x + / x] if < 0. 


This is quite a mouthful and we usually abbreviate it to 
e 
u(x + h) = u(x) + hu(x) +... + a u(x) + O(h"* 1), 


where O(h"* !) is read as (terms of ) order h"* +. Q(h"* !) means loosely that л" +! is the 
lowest power оЃл to appear in the expression it replaces. Note that О is not a function. 


More precisely, we define the notation as follows. We say that 
fih) = O(g(h)) as л 0, 
if J constants у > 0, K > 0 such that 
If (A) < Kigh) he[—n, 1). 
For example, if 
fü) = 3h, 
then, clearly we may say that f (h) = O(h?). In our case, we see (by Taylor’s Theorem) 


that the remainder in the nth Taylor approximation is O(h"*!) as h ~ 0 under the 
stated conditions, by choosing 


E Basi 
(n+ 1)! 


Taylor's Theorem isan important tool because it enables a function to be approximated 
in a given neighbourhood [x — h, x + h] of a point x by a polynomial in h, and 
polynomials are relatively easy to handle. In practice we shall use Taylor's Theorem 
without specifying that u satisfies the conditions. 


READ S: page 6, line 15 Finite-difference approximations to derivatives to page 8, the 
end of Chapter 1. А ' 


We have not asked you to read the first few pages of. dealing with elliptic equations 
(which are not discussed in this unit), but if you have time you might start reading 
from S: page 4 Parabolic and hyperbolic equations. d 


Notes 


ü) S: page 6, Equations (1.3), (1.4) and (1.5) 


The expression represented by O(h*) 
the left-hand side and the approxim 
error of the approximation since wi 


is the difference between the true value of 
ation. It is therefore called the local truncation 
€ have truncated the Taylor series to obtain it 
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The notation 
Е ` 
ах? |, 


Фаи(х) 


ах? 


теапѕ 


(ii) S: page 6, line —7 


The term involving the lowest power of Л in the error is the important one, 
because in practice h is much smaller than unity and we are interested in the 
behaviour as h approaches zero. 


(iii) S: page 8, Equations (1.8), (1.9) and (1.10) 


These equations are basic finite-difference approximations to the partial 


derivatives 
u e d ди 
ee ШИЙ ccs 
ex on дё 


evaluated at the point P with coordinates (x, t) = (ih, jk). (Note that the origin 

of our mesh of points is (хо, tọ) = (0, 0).) We shall make great use of these 

formulas throughout this unit in approximating the partial derivatives in a 
` differential equation. 


General Comment 


For a function of one variable, u: R —9 К, we have seen how forward differences 
may be expressed in terms of the forward-difference operator A,, whose effect on u 
is defined by 


Ayu zx кэ u(x + h) — u(x) xeR. 


With this notation, Equation (1.6) in S: page 7 can be rewritten concisely as 
, 1 
w(x) 5 {А„и(х)). 
1 


This is hardly surprising since we originally defined u(x) (Unit M 100 12, Differentiation 
I)as 


id 
lim = (A,u(x)]. 
n~o h 
However, it is pointed out in S: pages 6 and 7 that the central-difference approximation 
is more accurate. We therefore introduce the central-difference operator ó, on the 
set of functions R —> R, whose effect on u is given by 


Óju:x кэ u(x + dh) — u(x — $h). 


We may compose the central-difference operator with itself to obtain бу = à, ° Ôp, 
etc. Clearly 


ógu:x —» u(x + А) — 2u(x) + u(x — h), 
so that we may rewrite Equation (1.4) in S: póge 6 as 


l (i2 
iz (ógu(x)]. 


If we assume that the function и is tabulated at intervals of л, we encounter some 
difficulty in expressing Equation (1:5) in S: page 6 in terms of the central-difference 
operator б, since d,t(x) involves values of u at the untabulated points x + $h. We 
resolve this problem by introducing the averaging operator pp, whose effect on u is 
given by 


и (x) = 


pu ix > uix + S) + u(x — 3). 


We now have 
ju o! x i 3[бщ(х + Eh) + òpus 20) = Spi +) — ulx — A], 


and so Equation (1.5) becomes. with our notation, 


u(x) = 1 pó,utx)). 


EE 
h* 
Note the convention that the subscript on џ and the mapping product symbol are 


omitted when the meaning is clear from the context. 


When dealing with two independent variables x and t with, perhaps, different constant 
spacings h and k in the two directions, it is more convenient to suppress the subscripts 
hand k and write, for example, 


Ах 1) = u(x; + h, tj) — u(x;, t) 
Siu(x;, tj) = и(х + һи) — ux, t) + u(x; — h, 0). 
The subscript x here denotes differencing in the x-direction. 


With this notation as a useful shorthand, we can write Equations (1.8), (1.9) and (1.10) 
in S: page 8 as 


and 


respectively. 


These approximations have local errors, and we might wish occasionally to use 
more accurate formulas. The derivation of such formulas is quite lengthy and we 
Shall not deal with it here. However, we may, from time to time, quote results which 
you will have to accept on trust. (There is a good account of the theory involved in 
F. B. Hildebrand, Introduction to Numerical Analysis, McGraw-Hill, 1956.) 


10 
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5.2 EXPLICIT METHODS OF SOLUTION FOR PARABOLIC 
AND HYPERBOLIC EQUATIONS 


| READ S: the section entitled An explicit method of solution, pages /0. to: 17, omitting | 
‚ page 14, lines —18 to —7. 


We have not asked you to read the passage in S: pages 9 and 10 which shows that the 
one-dimensional heat equation for a finite rod can be expressed in terms of non- 
dimensional variables as 


2. 
би xe) E> 0 


Notes 
(1) S: page 10, Equation (2.4) 


We consider r to be a parameter whose value we can choose to make the formula 
relating the various values of и at the mesh points (or, as we shall call it, ‘the 
finite-difference scheme) as useful as possible. Another piece of jargon that we 
sometimes use is that we call the equation in S: page 10, line —5 the finite- 
difference replacement of the partial differential equation (2.3). The parameter r 
is known as the mesh ratio of the scheme. 


(ii) S: page 11, line 3 
It is usual to refer to a time level rather than a time row. The process of calculating 
values along one time level from the values on the previous time level(s) is 
called stepping forward in time where the step is of length k (the mesh length 
in the t-direction). 


(ili) S: page 11, Example 2.1 
Equation (2.4) on S: page 10 is valid only for points in the interior of the solution 
domain. We cannot make use of this formula to evaluate the solution at boundary 
points. However, we are given the boundary values along x — 0 and x — 1 and 
we are able to evaluate the solution at all points along a time level. 


(iv) S: page 12, line 8 
The choices for h and К (and hence r) are arbitrary at this stage. The only infor- 
mation we have about choosing values for /! and k comes from our discussion 
of Taylor's Theorem, from which we may conclude that "small" absolute values 
of h and k would be needed to give good results. 


Note also that there is an infinitude of choices for h and k which would yield the 
same value for r. For example, taking л = 1 and k = 7g also yields r = тр. 


(v) S: page 12, Fig. 2.2 
The idea of a molecule representing a finite-difference scheme is a useful one. 
Each atom represents the value of u at a single mesh point. The whole molecule 
gives a summary of the coefficients of a scheme and its form. (Note that a single 
atom at a future time level, i.e. above the "equals" sign, indicates an explicit 
scheme.) 


It also tells us how many time levels are involved, and how many initial and/or 
boundary conditions are required in order to use the scheme. This last point 
can be illustrated by considering an explicit scheme (on some level) of the follow- 
ing type, with suitable numbers in the atoms. 


1 
| 
! 
га 
га 
+ 
+ 
m 
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(vi) 


(Note that we have returned the "equals" sign to its familiar orientation. The 
: ‹ 
interpretation of the molecule remains thé same.) 


This, for example, could be produced by using. in place of the simple formula 


(You need not check this formula.) 


Because of the presence of two atoms to the right and left of i,j this scheme 
requires two conditions on each of the boundaries. Of course, we have these two 
conditions; namely, we are given 


(a) (0,)=0 120, 
say, and since this implies that 
ди 
= (0,) = 0 120, 
et ) 


we have, from the differential equation, 
Ou 
ax? 
Naturally the choice of finite-difference replacement cannot alter the number of 
boundary conditions required for the partial differential equation to be properly 
posed. 


S: page 13, line —6 
We shall ask you to obtain this solution in SAQ 21 of Unit 6, Fourier Series. 


(b) 00)-20  t20 


(vii) S: page 14, lines 1 to 8 


12 


The finite-difference scheme allows the effect of the discontinuity to propagate 
as an error in the solution. We can see this as follows. The points marked with 
a cross in the diagram are those points at which the numerical solution to the 
differential equation, using the finite-difference scheme (2.6) in S: page 12, 
depends on the value of u at (x, 0). We call this set of points the numerical domain 
of influence of the point (X, 0). 


d |] КИ 


tt (3,0) 


Now, for the differential equation which we 


en аге considering there is but one 
family of characteristics, t = constant, as we 


act | ^ have seen in Unit 2, Classification 
and Characteristics. Asa result, the discontinuity in ди/дх at (4, 0) does not cross 
any characteristic, and so does not propagate. 


On the other hand, we see that, for given h and k, the domain of influence of (s, 0) 
in the finite-difference scheme is bounded below by the lines 
k k 
t= Г — X) and t= -70 — X), 
along which the error propagates. (The speed of propagation is seen to be h/k.) 
For fixed h, as k approaches zero (and hence also r ~ 0) these two numerical 
characteristics both approach 


t=0, 
ie, the characteristic of the differential equation which passes through (x, 0). 


In fact, even if л is not kept fixed, but h and k both approach zero in such a way 
that г = k/h? remains constant, k/h = hr still approaches zero and the character- 
istics of the numerical scheme approach = 0. 


The error in our computation is generated by using a Taylor approximation to 
@?u/@x? at (5, 0)—a point where this second derivative does not exist! However, 
for our particular choice of h and k the scheme is stable with respect to this error. 
As the results in Table 2.4 show, the magnitude of the error decreases as t increases. 


The second paragraph of S: page /4 is referring to the numerical solutions of 
parabolic equations. As we shall not be investigating this effect any further 
there is no necessity to refer to Exercise 12 in S: Chapter 3 as suggested in line 8. 


(viii) S: page 15, Table 2.6 


(ix) 


The percentage error column of Table 2.6 indicates larger errors than in Table 2.3 
(S: page 13). This is the result of two effects. 


The first is that in Case 2 the rate of propagation of the error generated by the 
discontinuity is one fifth of that in Case 1, since now h/k = 20 rather than 100. 
We should therefore expect the results to be poorer until this effect has subsided. 
The second effect is that since k is much larger than before the local truncation 
error is larger. 


S: page 16, lines 1 to 4 

The values occurring in Table 2.7 which are less than 0 or greater than 1 obviously 
contradict the maximum principle for the heat equation (W: page 60, lines — 2 
and — 1). This gives us a warning that something is going wrong. 


S: page 16, line —3 

Smith describes a finite-difference replacement of a differential equation as 
valid if it is convergent, stable and compatible. These are topics covered in 
S: Chapter 3 and we shall discuss them in Unit 8. A short discussion of stability 
in relation to the explicit scheme of this section is presented in the following text. 


General Comment 


A possible reason for the failure of the numerical scheme in Case 3 arises from the 
fact that we are performing a step-by-step process with a recurrence relation (with 
two independent arguments і and j). We already know, from the work on recurrence 
relations with one argument (Unit M201 7), that such a computation can sometimes 
give rise to serious induced instability. ` 


In Table 2.7 (S: page 15) the instability does not arise through local rounding error 
because the entries in this table were computed with exact arithmetic as you can 
easily verify by hand. There is, however, another significant component of local 
error, present in all these computations. It is the truncation error which comes from 
the fact that Equation (2.4) (S: page 10) is not an exact replacement of the differential 
equation. Remember that each of the finite-difference formulas obtained in Section 
5.1 was subject to a local truncation error as we saw in note (i) of that section. The 
finite-difference equation (2.4), therefore, contains a local error which is a sum of 
terms O(h?) and O(k) appearing in Equations (1.8) and (1.10) of S: page 8. 


Formally, we define the local truncation error at a point of the finite-difference 
replacement of a differential equation as the ‘error it introduces in evaluating the 
solution at that point. i.c. the terms truncated in replacing derivatives by their Taylor 
approximations. For example. let U denote the true solution of 


subject to the stated subsidiary conditions. The local truncation error at i,j + | is 


I I 
Tu= rU = U;,) ue 


ur 2U;; + Ui. iJ. 


h? id 
If we can find the local truncation error at each point it will give us an indication of 
how well the true solution of the partial differential equation satisfies the finite- 
difference equation. In the present case we already know the true solution of the 
differential equation. Therefore, we have substituted the values ofthe true solution into 
the finite-difference equation for each of the three cases and have recorded the value 
of the local truncation error, T; ;, at selected mesh points. The results are shown in 
tables for the cases r = 0.1, 0.5 and 1.0 respectively. 


ioe S 


We can see from these results the different values of T; ; at any particular point in 
the three cases. For example, look at the figures corresponding to г = 0.02 (j = 20 
in Case 1, j = 4 in Case 2 and j = 2 in Case 3). We can account for these differences by 
the different sizes of К, in view of the O(k) terms in the local truncation error. 


"These results do not, by themselves, account for the drastic behaviour in Case 3. We 
can conclude that the effects of truncation errors made at each step are accumulating 
in the third case as the step-by-step process progresses. In the first two cases these 
errors must be diminishing in their effect. The only difference between the cases is 
the value of r and we conclude that it must be this which governs the stability or 
instability of the scheme. 


The large errors for j = 0, especially in Ts 9, are due to the discontinuity in ди/ёх at 
(3, 0). The fact that, in the first two cases, this error does not propagate through the 
solution illustrates dramatically how the effect of the local errors diminishes as t 
increases, for these cases. 


It is also interesting to point out that had Case 3 been computed on a digital computer 
where rounding errors are introduced, these too would have accumulated. However, 
rounding errors are often smaller than truncation errors, and in this unit we concen- 
trate on the latter. 


Note that in this section we have discussed only the parabolic equation 
ðu ди 
а ax 
There is an explicit scheme for the hyperbolic equation 
^u = Ou í 
Oe ax2 


Ox 


which may be developed along the lines discussed here. It has a similar property that 
its stability depends on its mesh ratio. This method is the subject of SAQ 1. 
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Case 1: = 01, К = 0.001, 7 = 0.1 
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і= 1 2 3 4 5 
х = 01 0.2 0.3 04 0.5 
(/ = 0) г = 0.000 2.23517Е — 05 2.08616E — 04 1.04308Е — 04 —7.88033E — 01 —3.13647E + 01 
1 0.001 |—4.47035E — 04 2.60770E — 04 2.50787E — 02 1.62046 — 3.67489 
2 0.002 2.23517E — 05 1.70246E — 03 1.77227E — 01 6.36525E — 01 — 1.63817 
3 0.003 1.49012E — 04 1.40928E — 02 245437E — 01 2.24337E — 01 —9.68307E — 01 
4 0.004 1.08778E — 03 346191E — 02 2.37897E — 01 5.45979Е — 02 —6.55293E — 01 
5 0.005 3.72529E — 03 5.29960E — 02 2.03758Е — 01 —2.01315bE — 02 —4.80771E — 01 
10 0.01 2.46763E — 02 6.81467E — 02 6.34789Е — 02 —7.27028E — 02 — 1.78874E — 01 
20 0.02 2.04742E — 02 242516Е — 02 —1.14739E — 03 —4.32655bE — 02 —6.47008E — 02 
Case 2: h = 0.1, k = 0.005, r = 0.5 
tel 2 3 4 5 
xl 0.2 0.3 04 0.5 
(= 0) г = 0.000 |—5.60284E — 04 — 3.05414Е — 02 —6.79231E — 01 —6.66516 8.08468 
1 0.005 |—6.26862bE — 02 — 3.49882Е — 01 —6.72478E — 01 3.26108E — 01 1.53010 
2 0.010 1.39344E — 01 2.92724E — 01 1.93632Е — 01 3.33082E — 01 6.92152E — 01 
3 0.015 |-1.26433E — 01 —1.75536E — 01 —3.60668E — 02 2.54286E — 01 4.13870E — 01 
4 0.020 |—8.60870Е — 02 —9.31948E — 02 2.13026E — 02 1.96314E — 01 2.82646E — OI 
5 0.025 |—4.97580E — 02 —4.05758E — 02 4.50074E — 02 1.57166E — 01 2.09284Е — 01 
10 0.050 1.50621Е — 02 3.53306E — 02 5.99980E — 02 8.13424E — 02 8.98718Е — 02 
Case 3: h = 0.1, k = 0.01, г = 1 
i=1 2 3 4 5 
х= 01 0.2 0.3 04 0.5 
(0 = 0) г = 0.00 |—3.89755Е — 02 —344896E — 01 —2.01015 — 7.98562 1.74325Е + 01 
1 0.01 |—3.42774Е—01 6.35606E — 01 —344312E — 01 7.68459E — 01 1.48830 
2 0.02 1.85341E — 01 1.88848E — 01 7.38025Е — 02 4.58241E — 01 6.44755E — 01 
3 003 |—4.47214E — 02 2.23517E — 03 1.42622E — 01 3.10755E — 01 3.85380E — 01 
4 0.04 1.68860E — 02 6.75470E — 02 1.53589E — 01 2.38597E — 01 2.74277E — 01 
5 0.05 3.94821E — 02 8.97855E — 02 1.48582Е — 01 1.98489E — 01 2.18272E — 01 
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5401 


ici 1 Р i al differences. fi 
Write down an explicit finite-diference scheme, using central differences. for thc 
hyperbolic equation ‘ 


Draw the molecule for your scheme. If the mesh length is / in the x-direction and 
k in the (-direction what is the mesh ratio (compare with note (i) of this section)? 
Write down the numerical domain of influence of some general point (X. г). Т > 0, 
for your scheme. What initial data does your scheme require? 


(Solution on p. 32.) 


In SAQ 2, which follows, parts (a) to (f) are designed to assist you in constructing a 
computer program which will yield the numerical solution to a particular case of the 
diffusion equation. If you are short of time then make use of the library program 
SEXP321, which has been developed as a solution to this problem, and tackle part (g) 
only. 


SAQ 2 


(а) Write a BASIC statement which can be used to evaluate шуу according to 
Equation (2.4) of S: page 10. Assume that иу, ti-1,j шуу; and r are already 
available in the computer. 


Р 


(b) Place the solution to (a) within a BASIC loop which will calculate u; ;,, for 
i= 2,3,..., N — 1 for some j, where N has already been read into the computer. 
Assume that и, (i = 1,2,..., N) are already available in the computer. 


(c) Incorporate the result of (b) into another BASIC loop which will calculate 
uj. forj = 1,2,..., M, where M is already available in the computer. Assume 
that ш; and uy; -.., M), and u(i = 2,3,..., N — 1) are already 
available in the computer. 


(d) In (c) what do the values 


(i) uanduy; j= 1,2,....M, 
{ шу F=1,2,...,N 


represent? 


(e) Why are zero values not used for the subscripts i and j in the previous parts of 
this question? 


(f) Write a complete BASIC program which will solve the initial-boundary value 
problem 


ди д?и zi б 
1 ae X«m is 


subject to the initial condition 

u(x, 0) = sinx О<хҳл 
and the boundary conditions 

u(0, t) = (лї) = 0 t20. 


Let N, the number of mesh points in the x-direction, and k, the mesh length 
in the t-direction, be read in as data, and arrange for r to be printed out. You 
may assume that the solution is not required for t > 20k. Theo 1 
program should show, at each mesh point: 


utput from your 
the analytical solution to the problem, which is given by 

u(x,t) = e^'sinx (x, 1) €[0, x] x Rg; 
the solution to the finite-difference scheme; 


the error, i.e. the difference between the two solutions. 


(в) 


PDE 52 


To avoid an inordinate amount of output show the analytical solution, numerical 
solution and error at only one value of x for each time step. Choose that value 
of x where you expect the error to be greatest. 


Run the program developed in part (f) with various values of N andk to investigate 
the nature of the numerical solution for various values of the mesh ratio, r, both 
greater and smaller than 1. Check your results by comparison with those 


obtained from the library program $EXP321. 


(Solution on p. 32.) 
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53 ANIMPLICIT METHOD OF SOLUTION 


5.31 The Crank-Nicolson Method 


READ S: the sections entitled Crank-Nicolson implicit method, pages 17 to 20, and 
A weighted average approximation, pages 23 and 24, omitting the final reference to 
Chapter 3. 


Notes 
(i) S: page 17. line 18 . | 
You may recall that we used this averaging device in Section 21.2.2 of Unit 


M201 21 10 obtain the trapezoidal-rule formula for ordinary differential equa- 
tions, and in Section 2.3 of Unit 2, Classification and Characteristics. The new 
formula has the property that the O(k) term of the truncation error for the explicit 
formula becomes О(К?) as we shall see in Section 5.5. 


(ii) S: page 18, lines —5 and —4 
Since k = rh? we must be careful not to choose too large a value for r, since the 
resulting large value for k would produce a large local truncation error. 


(iii) S: page 18, line —3 
The advantage in choosing r — 1 is that it simplifies the numerical formula and 
hence reduces the number of calculations required. This situation is of use in 
both hand and computer calculation since it reduces the total time required 
to obtain the solution. Note that we can still control the accuracy of the scheme 
with r = 1 by choosing an appropriate value for ћ (in which case К = А?). 


(iv) S: page 19, lines 5 to 9 
Note the matrix form of these equations: 


4-1 uy 0.4 
-1 4-1 Uz 0.8 
-1 4-1 из | = |12 
-1 4 -1)] ш 1.6 

-2 4| и; 1.6 |, 


which is of the form Ай = b where А is tridiagonal, that is, its nonzero elements 
appear on and adjacent to its main diagonal.* 


(v) S: page 23, lines —6 to —1 

This is an example ofa technique which is used widely in numerical mathematics. 
Having obtained a useful formula (in this case the Crank-Nicolson method) we 
try to generalize it to a formula which reduces to other well-known formulas 
for particular values of some parameter (in this case 0). Any analysis which needs 
to be applied to each of the particular formulas now needs only to be applied 
to the one general formula and the results for the particular formulas are obtained 
by substituting in the relevant value of the parameter. 


General Comment 


Once again we have restricted ourselves to a parabolic equation. Implicit methods 
can be found for hyperbolic equations. For example, for the equation 


* Such a matrix is called triple diagonal in Unit M201 30, Numerical Solution of Eigenvalue Problems. 
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an analogue to the weighted average approximation in the parabolic case is 


Ou 0n. (1 20) ёги, з + Оби 


ij IESU 
k? h? 


SAQ 3 
(a) Draw the molecule for the Crank-Nicolson method with mesh ratio r. 


(b) Write the Crank-Nicolson implicit method in terms of the central-difference 
operator 62. 


(Solution on p. 35.) 


SAQ 4 


Write a section of a BASIC program which will calculate the matrix of coefficients 
of the variables u; ;,,(i = 2,3,...,N — 1) generated by the Crank—Nicolson implicit 
method with general r. Let u, j and uy jj = 1, 2,..., M) be known boundary values 
and store the elements of the matrix in the array C. 


(Solution on p. 35.) 


5405 


Write down the matrix of coefficients for the Crank-Nicolson method applied to а 
problem involving six internal mesh points along each time level. Assume that the 
boundary values to j+; and из ууу are known. What special form does this matrix 
have? 


(Solution on p. 36.) 


The library program $CNH321 uses the Crank-Nicolson method to solve the equa- 
tion 


Ou — Qu 
Au a<x<b, Т> 0, 
e Ax? 


with the initial condition 

u(x, 0) = sin zx а<х<Ь 
and the boundary conditions 

и(ал) = 0 г> 0, 

u(b, t) = 0 г> 0. 


The program is written in BASIC. You can replace statements 160, 170 ог 180 by 
new DEF statements to specify alternative initial and boundary conditions for this 
problem. For example, the initial condition 


u(x.0) = x? + 6x +1 а<х<Ь 
requires that statement 160 be 

160 DEF FNI(X) 2 X* X + 6«X + 1. 
Input data consists of: 


G) the values a, b defining the boundaries (x = a and x = b) of the domain; 
Gi) the number of mesh points along each time-level; 

(iii) the value of the mesh ratio: 

(iv) the total number of steps to be taken in the r-direction. 


A listing of SCNH321 follows. The section between statements 400 and 500 solves a 
tridiagonal system of equations by the method to be explained in Section 5.32. 
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18 

24 

38 

40 

5B 

60 

76 

|20] 

9р 

100 
110 
128 
134 
140 
158 
155 
168 
178 
188 
198 
208 
210 
220 
236 
249 
256 
269 
270 
280 
290 
380 
316 
329 
ззй 
340 
350 
269 
378 
386 
396 
40И 
410 
42G 
435 
44g 
450 
2860 
476 
486 
490 
50g 
516 
520 
530 
54g 
550 
эби 
57Й 
580 
594 
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PRINT "THIS PROGRAM USES THE CRANK-NICOLSON E ес А 
PRINT "NUMERICAL SOLUTION ОЕ THE HEAT CONDUCTIO Boe ATENENT 
PRINT "THE INITIAL CONDITIONS ARE SPECIFIED BY A 
PRINT "AT LINE 168." 
PRIMI "BDUNDARY CONDITIONS ARE SPECIFIED BY DEF STATEMENTS TI 
PRINT "LINE 17Ø(FNA-B.V. AT X-A)AND LINE 180 (FNB-B.V. 
PRINT 
DIM A[22],8[22],c[28],D[28],G[20],U[ 28] , v [ 20] 
REMX**INPUT SEQUENCE*** А 
PRINT “INPUT END POINTS OF DOMAIN-X-A AND X=B"; 
INPUT A,B f 
PRINT "INPUT NUMBER OF MESH POINTS, N”; 
INPUT N 
H2 (B-A)/ (N-1) 
PRINT "INPUT PARAMETER R"; 
INPUT R 
DEF FNI(X)=SIN(3.14159*xX ) 
DEF FNA(T)=@ 
DEF ЕМВ(Т)=Й . 
PRINT "INPUT THE NUMBER OF TIME STEPS REQUIRED"; 
INPUT M 
PRINT 
REM***SET UP COEFFICIENTS OF LINEAR EQUATIONS*** 
FOR I-1 TO N 
A[I]-R 
B[I]=2+2*R 
C[I]=R 
NEXT I 
GOSUB 596 
PRINT ” SOLUTION” 
PRINT "ш ^ 
FOR JeB TO M 
GOSUB 650 
RINT "TIME- 
RINT "--------—- * 
F J=@ THEN 516 
MEX*CALCULATE R.H.S. OF LINEAR EQUATIONS*«x 
R I-2 TO N-1 
IJeR*U[I-1]«(2-2*n)*U[ T] -R*U[ T1] 
XT I 
M 
1 


22 отон о о 


E 
2 
E 
EM*SOLUTION OF TRIDIAGONA 


L SYSTEM* 


С[т]/(#&[11-А[т]%у[т-1] 
(ETI sal TI*el I-49] ) 7 (6 r]-A[r]*w[r-41]) 


REM*CALCULATE THE U(I)* 
FOR I-N-1 TO 2 STEP -4 
ULI]svu[I]*U[I-1]«G[I] 
NEXT I 
REM***OUTPUT OF RESULTS FOR ON Я 

Fu E TIME STEPxxx 
PRINT U[TI], 

EXT I 

PRINT 

PRINT " 

NEXT J 

STOP 
REM***SUBROUT IN 


E FOR INITIAL CONDITIONS ALONG Tagxxx 
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688 FOR I-2 TO N-1 

618 X=A+(I-1)#H 

c20  U[I]sFNI(X) 

630 NEXT I 

64g RETURN 

655p REM***SUBROUTINE FOR B.C.'S ALONG X=A AND X=B¥** 
668 T=J*H*H*R 

6798 U[1]=FNA(T) 

686 U[N]=FNB(T) 

698 RETURN 


764 END 

SAQ 6 

Use the library program $CNH321 to obtain the solution of 
Qu — hu 
u = x3 O0cx«l r»0 
Gt ex” 


subject to 

u(x, 0) = sin mx O<x<l 
and 

(0,0) = (1,0) = 0 120. 
at the point (0.5, 0.1). You should use a mesh length h = 0.1 and a mesh ratio r = 0.5. 


(Solution on p. 36.) 


5.3.2 The Solution of Tridiagonal Systems of Equations 


Although Gauss elimination with pivoting is a stable and natural method for solving 
sets of linear equations (Unit М201 8), the simple tridiagonal form of the system of 
equations generated by the Crank—Nicolson implicit method (see SAQ 5) allows a 
modification of the method which is simple to apply in practice. 


READ S: the section entitled Solution of the equations by Gauss's elimination method, | 
pages 20 to 23. 


General Comment 


The method of solution in the text is equivalent to Gauss elimination, but it turns out 
that no interchanges are necessary to ensure that the multipliers in the elimination 
are €1 in absolute value. As we showed in Unit M201 8. the elimination method is 
then stable. We should, however. prove the assertion that in our case the multipliers 
areall < 1 automatically, so that the pivots always lie on the diagonal of the coefficient 
matrix. 


In order to show that the method of solution of the equations in S: page 20 is stable 
we shall show two things. * 


(a) The multipliers in the process have absolute value <1, provided 
a; > 0, b; > 0,0; > Oand b; > aja, + ci. 
for each value of i. (We define cy = ay = 0.) 


(b) If in addition, b; > a; + с; for each value of i, then the back-substitution is stable. 
(We define a, = cy. , = 0.) 
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he sum of the absolute 


The final condition in (a) is that each diagonal element exceeds tl : ofthe: | 
їп (b) is similar with 


: ә p 
values of the other elements in the same column. The condition b) eee 
instead of "column". For symmetric matrices these two conditions are clearly 


“Tow 
equivalent. 

The matrix of coefficients for the Crank- Nicolson method satisfies all the conditions, 
since for each i 


a; = г, 

b; = 2 + 2r. 
and 

с=т. 


If you are short of time, you should omit the following proof. since it digresses from 
the main theme of the unit. 
Proof 


(a) From the first pair of equations in S: page 2/ we see that the multiplier m; used at 
this stage to eliminate u;.. is just а/о; _ 1, and from the formula for g; below Equation 
(2.12) we easily find that the multipliers m; = а/о; | satisfy the recurrence relation 


iet 


ты —L———— ї=2,3,...,М – 1, 
— cj jl 
with 
ay 
m, = =: 
b, 


We want to show that [m;, ;| < 1 for i = 1,2,..., N — 2, provided that 


a; > 0, b; > 0, c; > 0 and b; > aia, te 


Text 
We have 
ay : 
0<m, zm <1 since b,» 24,0. 
1 
Next 
аз , 
0 < m, =————— since a, > 0, b, > c, and m, < 1 
b, — сут Е T 
ay 3 
X ——— sinem < lande, > 0 
b-c, 
< eS I since Б 
г ee c) since >a [4 
(ау + с) — с, ? + 


and, by induction, it may be shown that 


О<т <1. 


(b) In Unit M201 8 it was asserted, but not proved, that if all |j < 1 then the back- 
substitution is also stable. We can in fact prove this easily for this 5р 


ion is als r ecial case, when 
the matrix is tridiagonal and no interchanges are needed in the elim 


ination. 
The required solution is computed from th К 1 
e backwards recurrence relati Д 
" A on 2 e 
21, line — 1) (S= page 


1 
t; = z% F Ciliga), 


where the S; and o; are given by the equations at the to 


i p of S: page 22. Nowi i 
7.32 of Unit M201 7 (on computations with recurr a cs 


ence relations) we saw that the 
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successive determination of u,.., %,—-2,-.., у from a recurrence relation of the form 
щ- = Piti + di. 


with u, specified, suffers neither from inherent nor from induced instability if |р < 1 
Гогі = п, п = 1....,2. 


Here p; = с; ,/o%;-,. and from the defining equations (ог g; we deduce the recurrence 
relation 


with 


Thus, using the given conditions, 
с, : 
0<р,=—<1 since b, > c, > 0. 
Next, 


since c; > 0, b; > a; and p, < 1 


since p < Landa, > 0 


< ——— = 1 since b, > az + С, 
(az + ¢2) — а 


and the general result, 0 < p; < 1, follows by induction. 
The method is therefore completely stable, and though the arithmetic is identical 
with that of Gaussian elimination, the recurrence relations in S: page 21, line —2 to 


page 22, line 3 express the method in a convenient and compact form for programming 
purposes. 


SAQ 7 


Calculate the solution to Example 2.2 of S: page 18 by the fully implicit backwards 
time-difference method (i.e. the weighted average approximation with 0 = 1) utilizing 
the recurrence relation version of the Gauss elimination method for the resulting 
tridiagonal system of equations. Take h = үр. and r = 1 and perform the calculations 
for the first time step only. 


(Solution on p. 36.) 


SAQ 8 


Write a BASIC program to solve a tridiagonal system of linear equations by the 
recurrence method. The input should include 


G) the number of equations to be solved, ; 
(ii) the nonzero coefficients, and 
(iii) the elements on the right-hand sides of the equations. 


(Solution on p. 37.) 
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SAQY i 
Solve the equations 
Эх, - ix, = l 
—3Хх| + xy - jx; = —2 
= х, + 2x, = dy, = -3 
-3 +20 =} 


by the recurrence method for tridiagonal systems using 


(i) exact arithmetic (use fractions). 
(ii) the library program STRI321 (see solution to SAQ 8). 


(Solution on p. 38.) 


SAQ 10 


Are the equations produced by the implicit scheme for hyperbolic equations (see 
General Comment in Section 5.3.1), for 0 0. soluble without instability by the 
method of Gauss elimination without pivoting? 


(Solution on p. 39.) 
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5.4 DERIVATIVE INITIAL AND BOUNDARY CONDITIONS 


5.4.1 Initial Conditions 


Consider the following initial value problem: 
Ou eu 
$7 md 0<х<1, t>0 


u(x, 0) = f(x) Oxxx«l 


ди 
= (x0) = (х) 0<х<1. 
д 
Suppose we use the explicit finite-difference scheme 
Unger = 203 + Uic. Mesa iij Щщ) 
k? h? 


tosolvetheinitial value problem. We can write the numerical scheme as 
Ны = 2(1 — ptg + pu; ъл, * i-i) — Hj, 


where the mesh ratio p? = k?/h?. We see immediately that to calculate values of u 
along t = (j + 1)k we need to know the values of u along t = jk and t = (j — Wk. 
The implication of this is that we shall need to specify initial values of u along both 
j = 0andj = 1. 


Wealready have u(x, 0) = f(x) givenalongt = Oso we use the other condition, namely 


Qu 
a 09 = g(x), 


to give us values along t = k. This derivative initial condition can be treated by the 
method of Section 5.1, that is, we replace it by a finite-difference formula. The choice 
of the forward-difference formula, 


ди (с.б) = (5) ыча 10, 
io 


yields the expression 


Ha — Uo = Кр 


where g; denotes g(xj). Now и; о = u(x;, 0) = f'(x;), by the first initial condition, and 
this allows us to write 


wa = + Ка  dP—0 d... 


where f; denotes f(x;). Thus we have obtained values along the time level г = k as 
well as along the initial line г = 0 and we can use the explicit finite-difference scheme 
to obtain the solution for t > К. 


5.4.2 Boundary Conditions 


Read S: the section entitled Derivative boundary conditions, pages 32 to 40, ignoring 
the references to Chapter 3. i 


Notes 


(i) S: page 33, line —6 
Beware! Here h is the positive physical constant introduced in 5: page 33, line 6. 
and по: the mesh spacing in the x-direction, which is now designated by àx. 
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S: page 35, lines — 10 and —9 
u is symmetric about x = 3 if and only if 
u(x, t) = u(1 х, t) O<x<l, 120, 
You may verify that if the function u satisfies the problem of Example 2.4 then 
so does the function 
(x,t ш = x.t) O<x<l 120. 

" . А : : i > 1 
Since the solution to the problem is unique, it must be symmetric about x = з. 
S; page 35, line —2 
Since the solution is symmetric we obtain цо; = toj, Wij = Ноу, Haj = Из, 
изу = Uy; and иу = Uj. and we need only solve for the six unknowns 
u; (i = 0, 1,...,5). The six equations required are those for t; у (i = 0, 1, 2, 3. 4) 
plus that for us j+; with ug j replaced by t, ;. 


54011 


Why is ди/дх represented more accurately by а central-difference formula than a 
forward-difference formula, as claimed in S: page 34, lines 3 and 4? 


(Solution on p. 40.) 


SAQ 12 


S: page 50, Exercise 5 


(Solution on p. 40.) 


You should omit SAQ 13 if you are short of time. 


SAQ 13 


S: page 51, Exercise 6 
If you are calculating by hand, compute и on just the first time level. 


(Solution on p. 40.) 
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55 ORDER OF LOCAL TRUNCATION ERROR 


In Section 5.2 we computed the local truncation error of a simple explicit method 
applied to the heat equation in one (space) dimension. That is, we found the extent 
to which the true solution of the differential equation did not satisfy the finite difference 
equation by computing 


1 1 
Ту = jun -U,j)- gi Ue 1j 7 20:3 + Uu 
where U; іѕ the value of the true solution U of the partial differential equation at the 
point (ih, jk). We have also implied that some formulas are locally more accurate 
than others, meaning that for the same / and k the local truncation errors T; ; are 
smaller in magnitude. 


We now distinguish between what we mean by the local error and the global error 
at a point. We view local errors in the following way. Suppose we know the true 
solution of the differential equation at all points up to and including those at time 
level j. Then if we were to use, for example, the simple explicit scheme to find an 
approximation u; j+ı to U; j+, we would calculate 


Unger = Uig + r(Uic i; — 205; + Uirigh 
where r = k/h?. 


By the definition of the local truncation error given above we see that 
Uj ja. Чужа = КТ 
measure the (local) error made in using the finite-difference 


In this sense the 7; 


a J 
scheme once only. 


In a complete step-by-step process we apply the finite-difference scheme many times; 
the local errors accumulate and it is this accumulation which produces the global 
error. Thus, at any stage, we can define the global error as the difference О; ; — tij 
between the computed finite-difference solution and the true solution of the differ- 
ential equation at a point. It should be clear from the unstable example in Section 5.2 
that simply reducing the local error does not ensure that the global error will be reduced ; 
the local truncation error, however, will provide a measure of accuracy when there 
is no accumulation of errors. 


We now show how to find an expression for the local truncation error T; ; for any 
finite-difference scheme applied to a differential equation. Obviously we cannot 
compute Т; in general because we will not always know the values of the true 
solution of the differential equation. What we can do is to use Taylor's theorem to 
give us an expression for Т;; involving the mesh spacings h and k and certain 


derivatives of the true solution of the partial differential equation. 


Again, in general, we will not know the values of the derivatives and at first sight we 
would appear to have gained no advantage. Fortunately, the expressions we obtain 
can be used to give a relative evaluation of various different schemes applied to a 
given differential equation. 


Before looking at a specific example note that we can write a finite-difference scheme 
in operator notation as 


ш-Ь=0 


where L is a finite-difference operator which replaces some differential operator and 
b is known at each mesh point and depends on initial and/or boundary conditions. 
The local truncation error, 7; у, is then the error which arises when the true solution 
of the differential equation is substituted on the left-hand side of the difference 
equation. Thus 
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To illustrate the complete method for finding T; ; let us consider the simple explicit 
scheme 


Шужа ОШ icq. — Ruj + Wiens, (1) 
k h? Ё 
applied to the partial differential equation 
eU — eu 
дг ax? 
in which the operator L = A,/k — 62/h? replaces 
8 0 
at ôx? 
We can write the local truncation error of scheme (1) as 
1 1 
T= gU -U- gu = 2U,,, + Uis, Q) 
Taylor's Theorem gives 
ôU; 1,,0?U,; 
Uji- Шу =k 3 i. 5k EZ 2 + Q(k) 


апа 


Огур 2U;;  Ui- 1; = (Uiig О) + (Ui-ig = Uig) 


Jj 1 0*0,; 
ij П] 4 ij 6 
ox? 12 дх* + 00 9 


under suitable conditions. Therefore, we can write Equation (2) as 
20, PU, kPU,; Y h? tU, 

jJ a ôx? ^2 at? 12 ax* 

Since U satisfies the differential equation we have 
QU; 220, 


ij 


ôt дх? 


J 4 O(k?) + O(h*). 


=0, 
and so 


“Од? 1 8x 
or, more simply, 
T; = O(k) + O(h’). 


4 + O(k?) + O(h*) (4) 


Sometimes in the literature the local error is calculated from the form of the difference 
equation given by 


Шужа = Mig + Tug — 20; t uus). 
That is, we would calculate 
Та = U 41 у Орр 2U,; + Шуу). 


It is easy to see that Tf, = KT, у-у. Either definition is acceptable and we shall use 
both in our numerical work. We have used the subscript i,j + 1 in our definition of 


T* because this is the local truncation error made in calculating u; ji 
гужа 


28 


PDE 5.5 


SAQ 14 


Show that the local truncation error of the Crank-Nicolson implicit method applied 
to the equation 


eU eU 
Ó ёх? 


is O(k?) + O(h?). 

(Solution on p. 41.) 

SAQ 15 

Find the local truncation error of the formula 


Шужа 205 ШШ, Mian 7 2unj T 


= 1.7 
k? h? 
when applied to the differential equation 
ги OU 
o! ôx? 


If the mesh ratio p? = k?/h? = 1 what value does the local truncation error take? 


(Solution on p. 42.) 


SAQ 16 


For what value of r — k/h? in the explicit scheme (1) would you obtain a much smaller 
local truncation error? (Use Equation (4) of this section as your starting point.) 


(Solution on p. 43.) 


a 
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5.6 SUMMARY 


s 


In this unit we have discussed and illustrated the following techniques: 


how partial derivatives may be approximated by finite-difference formulas 
using Taylor approximations; 

how to solve an initial-boundary value problem by replacing the differential 
equation by a finite-diflerence formula which is solved over a finite mesh 
of points in the xt-plane with corresponding initial and boundary conditions; 


the representation of a finite-difference scheme by a molecule; 


how to construct an explicit scheme for both the heat conduction and wave 
equations in one dimension; 


the Crank-Nicolson implicit method which leads to a tridiagonal system of 
equations; 


a direct method based on Gauss elimination, for the solution of a tridiagonal 
system of equations, which is conditionally stable and computationally 
efficient; 


how to include derivative initial or boundary conditions in the numerical 
scheme; 


the use of Taylor's Theorem to find the order of the local truncation error 
of a finite-difference scheme. 


We have also discussed induced instability in a numerical scheme due to the accumula- 
tion of local errors. For the explicit scheme this depends on the value of the mesh 
ratio. We have seen that a scheme determines the numerical domain of influence of 
a point and that the finite-difference replacement of a parabolic differential equation 
may have two families of numerical characteristics. 
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5.7. FURTHER SELF-ASSESSMENT QUESTIONS 


If you are short of time you may omit this section and return to it for revision purposes. 


SAQ 17 


S: page 46, Exercise 2 
You need not derive the analytical solution. 


(Solution on p. 43.) 


You may find the next SAQ somewhat tricky. 


SAQ 18 


Suggest an explicit finite-difference scheme for the diffusion equation in two dimen- 
sions, 

ди Ou ди 
6 Ox? ey? 
Use the same mesh spacing h in both the x- and y-directions. Draw the molecular 
diagram for the scheme. 


HINT: Draw a diagram (three-dimensional since и is a function of three variables) 
before devising the finite-difference scheme. 


As far as the actual computation is concerned, what added complications arise in 
this problem compared with the equation in one dimension, 


Ou @u Я 
à 0x7 


4 


(Solution on p. 43.) 


SAQ 19 


Describe a numerical experiment to solve the initial-value problem with which this 
unit began: 


e 


D Gern — C ur) + (sinxyu( t) = 0 0«x«l t»0 
дї дх? 


u(0, t) = (1,0) = 0 120, 
u(x, 0) = f(x) 0<х<1. 
(Solution on p. 44.) 
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58 SOLUTIONS TO SELF-ASSESSMENT QUESTIONS 


Solution to SAQ 1 


Using the central-difference formulas of Equations (1.8) 


and (1.9) in S: page 8 we obtain 


1 = Mij t з-д cag 7 uj + inj 


k? i 


gd 


as the finite-difference replacement of 


We can write the finite-difference scheme as 


Usa = 200 — pu, opus jj T i-i) 7 Meg 


where p? = k?jl? is the mesh ratio. The scheme. having only one unknown 4; 41+ 
is explicit. The molecule for the scheme is shown in the diagram. 


Bed 


ij-l 


For the numerical domain of influence we show all those points which depend for 
their value of u on the value at (X, Т) by crosses іп the diagram. 


NL. "A 
IN = х 
TN Zz 
NI IZ | 
LINE L. 


The numerical characteristics through (X, T) are 


—ї= p(x- X) 


and the numerical domain of influence is given by 


+, + фе —},—]+1,...,} — hjj = 0,1,2,...\. 


Since the scheme uses information on two levels (jand j — 1) we need to be able to 
specify values along г = 0 and t = k before the scheme can be used. For the Cauch 
problem these are available since both и and éu/ét are specified for t = 0. We st E 
see how this works in practice in Section 5.4. | Е 


Solution to SAQ 2 


(a) 


32 


50 LETU[LJ + I] = Uf, J] + R «(UL — 1. - 24 UL J) + U[L + 1, J) 


where U is a two-dimensional array which has been declared earlier in th 
program to be big enough to hold the solution to the differential equatio И F 
points in the solution domain. ation at all 


PDE SAQ 2 


(Generally speaking we would only require to store values of u along two time 
levels, j and j + 1, since once we have calculated и at all points along level 
j+ 1 the values along level j are no longer required. We could then replace 
ОП. J] by U[I, J + 1] and repeat the calculations to yield 


Uj j42 P= 1,2) 0.05.N.) 
(b 45 FORI=2TON-1 
50 LET U[L J + 1] = U[L J] + R *(U[I — 1, J] 2 2 + U[L, J) + U[I + 1. J]) 
55 NEXTI 
(c) 40 FORJ=1TOM 


45 FORI=2TON-1 
50 LETU[LJ-- t] = U[L J] + R«(U[L — 1, J] — 2s UÇ J] + ОП 41.4) 


55 NEXTI 
60 NEXTJ 

(d) (ì) u,;anduy;( = 1.2,..., M) are the boundary values. 
Gi) i — 1,2,..., № are the initial values. 


These values are assumed to be given so that they may be read in to the computer 
as data. ' 


(e) In BASIC we are not allowed to use the zero subscript. We have therefore 
assumed that the left-hand boundary is at i — 1 and that the initial data are 
given along j = 1. 


18 PRINT "THIS PROGRAM EVALUATES THE HEAT CONDUCTION EQUATION IN ONE" 
20 PRINT "SPACE VARIABLE USING THE SIMPLE EXPLICIT FINITE DIFFERENCE" 
30 PRINT "SCHEME" 

ай DIM UL295,211 

SØ PRINT "ТҮРЕ NUMBER OF POINTS IN X-DIRECTION AND SIZE OF TIME STEP" 
60 INPUT N,K 

70  Hz3-14159/CN-1) 

80  R-K/CHXHD 

90 М=20 

100 L=INTC(N+1)72) 

110  X-CL-10*H 

128 PRINT "VALUE OF R IS'5R 


130 PRINT 

140 PRINT " SOLUTION ALONG X="3X 

150 PRINT " J| ^-4-------------------2--2------- ы 
160 PRINT "Т1МЕ","ТНЕОН. VAL''5 "CALC. VAL « "", "ERROR" 

170 PRINT 


180 FOR I=2 TO N-1 

190  ULI»12J-SINCCI-10*H) 

900 NEXT I 

210 PRINT f5SINCCL-12*H25»ULL51150 
220 PRINT 

230 FOR J-1 TO M 

240 ULt1»Jl1-UEN5»J1-0 

250 NEXT J 

260 FOR J-1 TO M $ 

270 T=J*K 

280 FOR I=2 TO N-1 

2990  ULI,J*11-ULISJItR*CULIT15 J1-28*ULI» JI+ULI-12J]) 
300 NEXT I 

310 A=EXPC-T)*SINCX) 

320  C-ULL5,J-*11 

336 PRINT Т›А›,С›,А-С 


340 PRINT 
350 NEXT J 
360 END 
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The preceding listing is the program SEXP221. 
(g) The following results are a typical set. obtained with N = 5, k = 0.1. 


RUN 
EXP321 


THIS PROGRAM EVALUATES THE HEAT CONDUCTION EQUATION IN as 
SPACE VARIABLE USING THE SIMPLE EXPLICIT FINITE DIFFERENC 


SCHEME p 
TYPE NUMBER OF POINTS IN X-DIRECTION AND SIZE OF TIME STEP?55 0.1 


VALUE OF R IS -162114 


SOLUTION ALONG X= 145708 


TIME TH EOR* VAL САС,» VAL» ERROR 

@ le le Ø 

ol ° 994838 905036 -1.98245Е-04 
*8 +818731 81999 -3. 58939E-04 
23 » 740818 ° 741306 -4.867328E-04 
4 «67032 + 670908 -5.87940E-04 
.5 ‚ 606531 «607196 -6. 65188E-04 
„6 «548812 549534 -7+ 22403 E-G4 
7 «496585 ‚497348 -7.62463E-04 
8 © 449329 ‚450118 -1-88450E-04 
„9 40657 ©407373 -8.02755Е-04 
1 •367879 •368 687 -8.07285Е-04 
1.1 332871 «333675 -8.03471Е-04 
1.2 +301194 ©3019387 -7.93219 E-04 
1.3 „2725382 .273309 -7+ 77602E-04 
1.4 | «246597 ` +047355 -7+5761ЗЕ-04 
1.5 ^.22313 «223865 -7.34776Е-04 
1.6 .201897 202606 -7-O9146E-G4 
1-7 182684 +183365 


-6.81877E-04 


1.8 +165299 +165952 - 6+ 5338 6E-Q4 
1.9 ‚149569 .15g193 -6.84890E-04 
2 -135335 -13593 -5-94527E-04 
DONE 


34 


PDESSAQ3 4 


Solution to SAQ 3 


(a) The Cra 


nk-Nicolson method. given in S: page 17. is 


Sige (242и рр rus gua mr VO 20) ц jt rus, у. 


Its molecule is 


242r 
i-Lje-l jt КАНЕ! 
2-2, 
© О = 
i-dj ij + 


(b) We have 


and 


<2 
душу = uua; T ij + iag 


52 
ORM je m Моята lijer t Mena jen: 


Therefore, we can write the Crank- Nicolson method as 


or 


“52 ee A 
Quy pe — Гб ржа = ORM; + 2и, 


(2 — r2 jay = (2 + гб). 


Solution to SAQ 4 
Let the two-dimensional array U[LJ] (I = 1,2,..., N; J = 1.2,..., M) hold the 


values of the 
The coefficie: 


solution of the finite-difference scheme at the mesh points (x, t) = (ih, jk). 
nts of the unknowns at the (j + 1)th time level are to be stored in the 


(N — 2) x (N — 2) array C. 


We assume that the arrays C and U have already been declared and that the solution 
to the differential equation is known at all time levels up to and including г = jk. 

50 FORL-2TON-3 

60 LETC[LL —]] = -R 

61 ТЕТ СШ) = 2 + 2* К 

62 LETC[L L 41] = -R 

70 NEXTL 

80 LETC[i,1] =2+2+#R 

81 1ЕТС[1,2] = -R 

85 LETC[N-2,N-3]- -R 

86 ЕТ С[М – 2, М – 2] =2 + 2*R 

The boundary values, which are known, are assumed to occupy (1. J] апа U[N, J] 
(J = 1.2,.... М). Asa result there are N — 2 unknowns to be found at each time level. 
(Note that the piece of program in this solution is inefficient. In practice we would 


first calculat 


е2 + 2r and —r and use stored values for the LET assignments.) 
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Solution to SAQ 5 
The Crank- Nicolson formula is i 
Т =r 2 = 2r)u;; + ery 
rier ger t Oo 2r jua Z isije = Doa T(Q2-2ru; + Г; 


i i і = О and і = 7 = he 
For six internal mesh points let i — 0, 1.2,....7 where i =0 and i =7 are th 


boundary points. The relevant equations are 
HM jer t (2 + 2r jua maja = b, 
arty jua cb (2 + 20) jer зна = b; 
rita jar + (2 + 2a jii 7 Maja = b, 
—rus jay c (2+ Anus gea 7 Ms je = b, 


srg jay Qo 2rus ua — Moje = bs 


t+ + + + + 


rts jay + (2  2r)ug jer — 05а = by 


where Бу, ba. b3, by, bs and b, are obtained from the known values of u along the 
jth time level. Hence, in matrix form, 


242r —r LINES b, + моъи 
-r 242r -r Us j+ b, 
-r 242r -r Us jer b, 
=r 242r -r Ug jet Е b, 
=f 24+2r =r || bs 

=r 24 2r] [ue jar be + rus jad 


where no j+1 and ну; аге known boundary values. 


The matrix is tridiagonal and also symmetric. 


Solution to SAQ 6 


At (x. 1) = (0.5, 0.1) the solution of the problem using the Crank- Nicolson implicit 
scheme with r = 0.5 and Л = 0.1 is 0.3756, with a percentage error of 0.8. This is 
better than the explicit method with r — 0.5 but not so good as the explicit method 
with r = 0.1. (Further results are presented in S: pages 49 and 50.) 


Solution to SAQ 7 
Putting 0 = 1 in the weighted average approximation gives the scheme 
Wei T Hj = (Ura y jui = Riti jsa т Uj- asl) 
which can be written as 
"алуа t Qr duisi Mj m y 
The values of u for the first time step satisfy 
3u, — uy, = 02, 
—u, t 3u,— uy = 04, 
=u + 3u4 — uy = 0.6, 
—us + 3u — us = 08, 
= 2u, + 3u, = 1.0, 
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where и; denotes и, у. To obtain the last equation we have taken advantage of the 
symmetry of the problem. Using the notation of 5: page 20 we have 

d;—78,-d04,-—l.a5— 2, 
b = b, = Бу = Ба = b; = 3, 
E 
d, = 0.2, d; = 0.4, d, = 0.6, d, = 0.8, ds 


1.0. 


We now define g; and S; recursively, for i = 1, 2, 3. 
of S: page 22. 


5, using the equations at the top 


1 
t; = —(S; + cuui). 


we obtain 
us = +24 = 0.821 
ug = 39 = 0.732, 
uy = 28 = 0.574, 
ua = 42 = 0.390, 
и, dH = 0.197. 


Solution to SAQ 8 


The following is a print-out of the library program STR1321 which uses the recurrence 
method for the solution of a tridiagonal system of linear equations. 


PRINT "PROGRAM TO SOLVE A SYSTEM OF TRIDIAGONAL LINEAR EQUATIONS" 
PRINT 
REM 


16 
29 


Зӣ X*** THE INPUT STAGE  *** 


PRINT 
INPUT 
PRINT 
PRINT 
PRINT 
INPUT 
FOR 
PRIN 
INPU 
NEXT 
PRIN 
INPU 
PRIN 
MAT 
REM 
5[1] 


"ТҮРЕ THE NUMBER OF EQUATIONS(<=19)"; 
N 
"TYPE 3 SIGNIFICANT ELEMENTS OF TRIDIAGONALS FOR EACH ROW" 
= (N.B. ONLY 2 ELEMENTS FOR FIRST AND LAST Rows)” 
"FIRST ROW"; 
8[ 11, C[ 1] 
I=2 TO N-1 
T "ROW";I; 
T A[r],B[I],.C[I] 
I 
T "LAST ROW "; 
T A[N]. BIN] 
T "TYPE R.H.S8. 
INPUT DN] 
+ CALCULATE THE L(I) AND S(I) 
=р[ 1] 


CONSTANTS" 


++ 
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208 
218 
220 
230 
240 
25g 
260 
276 
288 
29g 
388 
318 
32g 
ззй 
348 
358 
36g 
376 
388 


L[11-8E 1) 

FOR I=2 ТО М 

Z-A[I]/L[I-1] à 
L[r]eB[r]-c[Ii-1]*z 
5[1]=0[1]-511-1]*2 

NEXT I 

REM *** CALCULATE THE U(I) *** 
U[NJsS[N]/L[N] 

FOR I-N-1 TO 1 STEP -1 
u[rJe(s[r]-c[r]*ufr«11)/L[I] 
NEXT I 

REM 3** OUTPUT STAGE *** 
PRINT 


FOR I=1 TO N 
PRINT U[I] 
NEXT I 

END 
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Note that the arrays store the actual coefficients, so that if a; and c; are defined as in 
S: page 20 then A[I] and C{I] store the values —a; and — су. 


Solution 


to SAQ 9 


(i) Using the notation of S: page 20 we obtain 


7 
d, = -1,d, = —2,d4 = —3, d4 = 
Непсе 
a = у = 2, 
аз 313 7 
P UTR UNA 
M T 222 % 
=. 8, 29 383 4 
as b; 2227357557 
а, 373 95 
a4 = by — Że =2 + 2:2: 2, 
ата 2427 16 
and 
S,;=d,=-1 
$,2d 49352. 3 3.1, H 
9, 22 4 
dg aegri ilu 54 
a3 2 7.4 
7 37 54 95 
S,=d + #$ mq улы a но 
65$ 973* 3377773 
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Therefore 
жы Ж у С 
© a 4 95 
1 И 7 54 3 
х = 5 (Ss жох) ci Fada a3 


(ii) The output from the program $TRI321 for this problem follows. 


RUN 
TRI321 


PROGRAM TO SOLVE A SYSTEM OF TRIDIAGONAL LINEAR EQUATIONS 


TYPE THE NUMBER OF EQUATIONS («217 )?4 

TYPE З SIGNIFICANT ELEMENTS OF TAIOIAGONALS FOR EACH ROW 
(N.B. ONLY 2 ELEMENTS FOR FIRST AND LAST ROWS) 

FIRST R0W?2,-1.5 

ROW 2 ?-1,5,2,-1.5 

Row 3 2-4 552.7125 

LAST ROW ?-1.5,2 

TYPE R.H.S. CONSTANTS 

?-1,-8,-3,3.5 


SOLUTION 


DONE 


Solution to SAQ 10 

The implicit method for the hyperbolic equation is 

бф; = p(0ó2u уа (1 — 20)ó2u; j + 052и, ;- 1) 

where p? — k?/h?. Assuming that the values along time levels j and j — 1 are known 
this gives rise to the equation 


2 
Ujei—p Olti jsi — Qu; ja + Usage) di, 


where d; depends only on known values. If we assume that the boundary values 
Hg, jas and иу j+ are also known, then, combining the equations fori = 1,.... N- 
we obtain ] E 
250-1  -p0 uy d, + p?0ug 
—р?0 2р0+1 -p0 и» 4, 
-p0 2р0 +1 —p?0 uy | =|d3 
=р?0 = 2p?0 + 1] |у. dy, + p?Ouy 


where we now write u; for j+,- It is clear that, in the notation of S: page 20, 
a; > 0, b; > 0, c; > 0, b > aj + c; and b; > ау + сл: Therefore, the recurrence 
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method for solving tridiagonal systems (which is equivalent to Gauss elimination 
able when applied to this system. 


without pivoting) i 


Solution to SAQ 11 


We saw in S: pages 6 and 7 that the central-difference approximation to a first 
derivative has an error of order А? whereas the error in the forward-difference 
approximation is of order h. Since h is less than unity, the conclusion follows. 


Solution to SAQ 12 
We use the explicit scheme 

шыжа = Hoy + Tau ij + Mia) 
(Equation (2.4) in S: page 10). 


(а) For the internal mesh points (i = 1.2..... N — D) the explicit scheme may be 
used directly. At the boundary x = 0 the condition 


= h = r) 


ё 
ёх 
may be approximated by the central-dilfcrence formula 


і 
3j Und — uoa) m e; = Гу). 


Elimination of v_, ; between this equation and 
Uo jer = Uoj + (Ирр oj + и.) 
yields 
Hos m (1 = 2r(0 + hi Òx) Hoj + wuj + 2rh vx. 
Similar treatment of the boundary x = 1 yiclds 


Uy jay = 2ruy aj t d = 2r(0 + hx) yuy у + Bhyvydx. 


(b) At the boundary x = 0 the condition 


ĉu 
ex 


= (и = vy) 
may be approximated by the forward-difference formula 


x a = uo.) = (ноз 04). 


The scheme yields 


Usa = Uj o, — 2u, j + uy) 


{1 — 2r + r + WyOx)}uy + rus + rhydx/(t + Nox), 
and the forward-diflerence formula may be applied at level j + 1 to give 
Ho. ji = (jay + hv dx)/ + Bx). 


The equations for uy- 1j+1 and uy j+, are obtained similarly. 


Solution to SAQ 13 


Modelling the problem is straightforward. The solution to the problem 


appears i 
S: pages 51 to 53. You need not verify the analytical solution. ass 


The function erfe in the solution arises as follows. Many functions have indefinite 


integrals which cannot be expressed by the use of a finite number of elementar 
e jals «i 3 Y " i 
functions (polynomials, sin, cos, exp, log). Many of these Integrals occur so even 
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that they have been given names and have been tabulated as new functions. One 
such function is the (so called) error function defined by 


The complementary error function is defined by 
De [9 oos 
erfe(x) = 5 ede xeR, 
n x 


from which it can be shown that 
erfe(x) = 1 — erf(x) xeR. 
Values of the functions erf and erfc are tabulated in, for example. M. Abramowitz 
and I. Stegun (eds.), Handbook of Mathematical Functions (Dover, 1965). 
Solution to SAQ 14 


The Crank-Nicolson formula is given on S: page 17, and the local truncation error is 


1 | ged 2 

Ty = guises = Ug = эз 8:01 + буй, )). 
where U satisfies 

eU _ eu 

e OX 
Now, by Taylor's Theorem. 

OU; 020; 3 
Uz jar = Uij + k P + 2 à d + O(k?) 


under suitable conditions. 
As in Equation (3) of Section 5.5 we have 


,8U; P OU, 


J + Q(h5), 
ae tpu ae 5000 


20, = 


апа 
a2 4 a4 
28 Шы+‹ 1 0,1 


ax? 12 óx* 


BU ja =h + O(h*). 


Applying Taylor's Theorem to ¢?U/@x? and @*U/éx* we obtain 


d hj 
ОУ zit 
ox сх? 
and 
OU ey CFU, 
je ij 
= = — 3 + 000) 
ex* ex* 
Thus 


: ae = c | - 2 Шы + O(k?) + O(h*) + О(КА?). 
Since U satisfies the differential equation, 
et ex? 
and the local truncation error is 
Tj = = шы + O(k?) + O(h*) + O(kh?) 


O(k?) + O(h’). 


Ш 
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Solution to SAQ 15 


Expansion in Taylor scrics yields 


7 kt бё, 
Uy jay = 203 + Urai = s Di 
and 
he OU; те 
Оны 360 ox 
Therefore, 
HS CUL І 20% p EU 
HC ae 12 art üx* 
1 O°U;; Q*U,; 
— [k*—— pm... 
360 | are дх® 
The differential equation gives 
д2 az 
сы eli 
et Üx^ 


Lf UU) 1 рги, aU 
a „2 Uij _ p28 Ui i “Шы _, j 
Т, l ae a) + 560 ре СЛ aet 
and the local truncation error of the scheme is O(k?) + O(h?). 
If p = 1 then k = hand 


т [2 x Lu ы ы) 


у= fi) p Em 
"A 421 er ox 360 | д дх° 

Now, assuming that U has partial derivatives of all orders, 

eu 2o s 


at ay? 


ac] using the differential cquation 
Ox" 


д? [o?U А n ; 
|= since all derivatives exist 


A 4 using the differential equation 


We may show, by induction, that 


ĉu е2" 
Ge = уат п = 1,2,... 


whence Т; ; = 0. We have shown that the finite-difference scheme 
Шужа = Mpeg, T Mica m Wigs, 


is an exact finite-difference representation of the wave equation. 
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Solution to SAQ 16 


The local truncation error, given by Equation (4) of Section 5.5, is 


kPU, Meu, 
Te ace eed ij 2 4 
е5 3p - Fa a ОК) + О) 


Provided that the third derivatives of U all exist and are continuous at i,j we have, 
using the differential equation, 
ёи mu eu eu 
e — Q8! O00 дх* 


at ij, 


and we can reduce the error to O(k?) + O(h*) by choosing r = k/h? = 1. 


(This is not very useful, however, since if h is small then k = $h? is very small and the 
resulting volume of computation is quite substantial.) 


Solution to SAQ 17 


See S: pages 48 and 49. i 

When it is not possible to remove a discontinuity by finding a transformation as 
suggested in S, it may be feasible to use very small values of the mesh spacings for 
the first few steps to improve the accuracy until the effect of the discontinuity has been 
sufficiently reduced to enable larger, more economic values to be used. 


Solution to SAQ 18 


7—77 с = (п + Ik 


Ры ЕТ” АР” АЕ” E ЭЙР” 
L ELL LLY’ 
149787 49 A AP P A AVIA 


LZ 
PP AP ARP APA E E A 
22212712 271207 


vw 


Z 
КПА V IL эй 

КГ 2 2 

CL La 52,2 


3 ғ 1 = пк 


The diagram shows two time levels, t = nk and t = (n + 1)k, where k is the mesh 
spacing in the t-direction. We replace the space derivatives by central differences 
and the time derivative by a forward difference. Using the notation 


Yijun 


to represent u(ih, jh, nk), we obtain 


Wise] — Hia _ Шыужыл T Mija Ж Hijo tn a Ura jn — Mijn + Ho in, 


k ПЫ h? 


which can be written as 


Un jms = Mija E Tei T Uij-in + Меваи T Шаа lli jn) 


where r = k/h? is the mesh ratio. Since there is only one unknown at the (n + 1)th 
time level the scheme is an explicit one. 
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The molecule of this scheme is shown in the diagram. 


i+ hin 


There is no difference, as far as the computations are concerned, between this case 
and the problem in one space dimension. There is, however, a greater number of 
mesh points to be taken into account. Note also that, provided we take a square 
mesh in the space dimensions, the mesh ratio is the same in both cases. 


If we were to consider implicit methods then we should have to solve (N — 1} 
equations in (N — 1)? unknowns at each time level (assuming a square mesh with 
N — 1 internal mesh points in both the x- and y-directions). The problem, for higher 
dimensions, is finding computationally faster methods for solving the much greater 
number of simultaneous linear equations which arise in the case of implicit methods. 


The Crank-Nicolson implicit method gives 
"оҳ 2 
Higanti ja = 5 ID + бушана + бщ jn + ju, a]. 
Solution to SAQ 19 
We can write down the explicit numerical scheme 


Ujar Чы Mang 20р nay Р 
> = 3 Я = — sin(ihjk)u; j 
n 1 (м, ; 


or the Crank-Nicolson implicit scheme 


2 h? k? 


Ujei — Mig I fuisi;s1— Qui jaa + Щ-1у+1 Q Hii 7 2j + ej 
k 2 Б 


= {зілі + 1)k)u; jay + sin(ihjk)u, у}, 


both of which can be solved by the techniques covered in this unit. The initial condition 
becomes 


uno = J (ih). 


The only problem is the possibility of a discontinuity at (0,0) and (1,0) if f(0) or 
/(1) is nonzero. The effect of such a discontinuity can be reduced by taking small 
mesh spacings for the first time steps. 


Unit 6 Fourier Series 


Contents 


6.0 


6.1 


6.2 


6.2.0 
6.2.1 
6.2.2 
6.2.3 


6.3 


6.3.0 
6.3.1 
6.3.2 
6.3.3 

6.4 


6.4.0 
6.4.1 
6.4.2 
6.5 
6.6 


6.7 


Set Books 
Conventions 
Bibliography 


Introduction 
Separation of Variables 


Convergence of Fourier Series 


Introduction 

Uniform and Mean Approximation 
Parseval's Equation 

Dini's Test 


Generalized Trigonometric Series 


Introduction 

Complex Fourier Series 
Sine'and Cosine Series 
Change of Scale 
Applications 


Introduction 

The Heat Equation 
Laplace’s Equation 
Summary 


Further Self-Assessment Questions 


Solutions to Self-Assessment Questions 


Set Books 
G. D. Smith, Numerical Solution of Partial Differential Equatious (Oxford, 1971). 
Н. F. Weinberger. A First Course in Partial Differential Equations (Blaisdell, 1965). 


It is essential to have these books: the course is based on them and will not make 
sense without them. They are referred to in the text as Sand W respectively. 


Unit 6 is based on W: Chapter IV, Sections 14 to 16, 18, 20 to 23. 


Conventions 
Before working through this text make sure you have read A Guide to the Course: 


Partial Differential Equations of Applied Mathematics. References to Open University 
courses in mathematics take the form: 


Unit M 100 13, Integration І for the Mathematics Foundation Course, 
Unit M201 23, The Wave Equation for the Linear Mathematics Course. 


Bibliography 
E. M. Horsburgh (Ed.), Napier Tercentenary Celebration (Royal Society, 1914). 
This book was produced on the tercentenary of Napier's birth. Chapter 6 gives a good 


description of various harmonic analysers, describing both the mechanics and the 
mathematics. 


J. Harvey, ‘A Harmonic Analyser’ in Proceedings of the Physical Society, 1930, Vol. 42, 
pp. 245-49. 


This is a paper describing a mechanical device for analysing a function into its 
harmonics. 


R. Furth and P. W. Pringle, ‘A New Photo-Electric Method for Fourier Synthesis 
and Analysis’ in Philosophical Magazine, 1944, Vol. 35, pp. 643-56. 


This paper described an electronic harmonic analyser. 


M. Spivak, Calculus (Benjamin, 1967). 


This is the set book for the Analysis course, M231. 


J. Fourier, The Analytical Theory of Heat (Dover, 1955). 


This translation of Fourier’s original work is primarily of historical interest. 


6.0 INTRODUCTION 


Our main aim in this unit is to revise one of the most common analytical methods of 
solving partial differential equations. The basic ideas have been covered before, in 
particular in Unit M201 22, Fourier Series and Unit M201 32, The Heat Conduction 
Equation. although a few of the concepts are new. 


The first part of the unit (Section 6.1) deals with the method of separation of variables. 
Separability is an extremely important property of many linear partial differential 
equations, as it allows them to be transformed to systems of ordinary differential 
equations. In general, analytical solutions are available for separable problems only. 


The ordinary differential equations that are obtained from separating the variables 
lead us to look for solutions which can be expressed as Fourier series. 


Thus the next part of this unit comprises a revision of Fourier series, which were 
introduced in Unit M201 22, and a short study of their convergence properties; to 
do this it is necessary to introduce a few notions of real analysis, in particular the idea 
of uniform convergence. ` 


Finally the ideas described in the earlier parts of the шій are used in a few of the 
more common problems in physics. 


You will notice the large number of SAQs in this unit. We have set these deliberately, 
for a considerable amount of the material contained here is revision. We expect you 
to concentrate on those SAQs which are in the areas where you feel you need most 
practice, and not to bother with SAQs which test material that you are confident 
about (although there is no harm in glancing through the solutions even in this case). 
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61 SEPARATION OF VARIABLES 


READ W: Section 14, pages 63 to 69. 


Notes 


(i) 


(ii) 


(iii) 


W: page 68. line — 17 | : 
The idea of a "closed form" is rather vague, and does not admit a rigorous 
definition. However. it may be taken to mean an expression involving the 
elementary functions which can be evaluated using a finite number of arithmetic 
operations. 

W: page 68. line — 15 | 
We shall be studying the relevant part of Chapter VII in Unit 13, Sturm-Liouville 
Theory. 

W: page 69, lines 8 and 9 

This condition amounts to 


ди 
«(у)ц(0, y) + [1 — «(] 2:00, у) = 0 


with «(y) = 1 for some values of v and о(у) = 0 for the remainder. Thus the 
coefficients are not independent of y. 


General Comment 


We have shown that if an equation is separable a solution can be written as а product 
of functions of one variable. This is obviously of computational convenience, and 
so it is clearly important to know when an equation will admit a separable solution. 
It is interesting to note that there is no simple rule for determining whether this 
is so, and that the best that can be done is to try different variables in turn. 


5401 


Assuming that the equation 


Oy öf ду 
abe 
Ot^ ax \ 6x 


has a solution of the form y = X (x) T(t), determine the equations satisfied by X and T. 


(Solution on p. 27.) 


540 2 


You obtained the two-dimensional Laplace operator in polar coordinates in SAQ 5 
of Unit 3, Elliptic and Parabolic Equations. Use it їо. па a solution of 


DV O<r<l, 


(where r is the radial coordinate) which 


(i) 
Gi) 


possesses radial symmetry; 


is bounded; 


(11) has the value 1 on the unit circle r = 1. 


HINT: Look for a solution of the form $ = r*. 


(Solution on p. 27.) 
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6.2 CONVERGENCE OF FOURIER SERIES 
6.2.0 Introduction 


Separation of variables leads us to believe that in many cases the solutions to problems 
involving partial differential equations can be expressed as infinite series. This type 
of solution first appeared in analysis in connection with the investigations by Daniel 
Bernoulli (1700-1785) of vibrating cords. He showed in 1748 that a formal solution 
to the problem of determining the motion of the cord is 

= .OHRX nct 

Y= У b,sin 7 соб 170, 

п=з I l 
the ends of the cord being fixed at x = O and l; and he asserted that this was the most 
general solution to the problem. 


A criticism of Bernoulli's theory was published immediately afterwards by Euler, 
who pointed out that a consequence of the theory was that any function f with 
f (0) = f() = 0 could be expressed in the form 


Е . ATX 

Ло) = Y b, sin TE 
net l 

To Euler this was absurd since the series is odd and periodic whilst f (x) need not be. 
This disagreement led to enquiries into what functions really were and the subsequent 
investigations touched upon the very foundations of the idea. Eventually this pro- 
duced the formulation with which you are familiar from the Mathematics Foundation 
Course. 


In the eighteenth century, f was called a function only if f (x) could be represented by 
a simple analytic expression such as a polynomial or a power series. But a mapping 
with an arbitrary graph, e.g. a polygonal line, was not accepted as a function. The 
existence of Fourier series, however, showed that such graphs could define functions. 
It was a long time before this matter was clarified, and the resulting studies were 
instrumental in the understanding of the foundations of analysis. For example, 
in 1861 Weierstrass used Fourier series to give an example of a function continuous 
throughout its domain but without a derivative anywhere.* 


Fourier series take their name from Jean Baptiste Joseph Fourier (1768—1830) who 
based on them his mathematical theory of the conduction of heat presented to the 
French Academy in 1807. Here he developed the theory of Fourier series and their 
application to boundary value problems in partial differential equations. He also 
enounced the proposition that an arbitrary function given graphically by means of a 
curve, which may be broken by ordinary discontinuities, is capable of representation 
by means of a trigonometric series. This theorem is said to have been received by 
Lagrange with astonishment and incredulity. 


After this work it was generally agreed that Fourier series were respectable. However, 
there still remained a few heretics, since Fourier did not provide a proof that the 
series actually converges at each point to the value of the function concerned. This 
point was not resolved until almost a century had elapsed. 


The importance attached to Fourier series in applied mathematics is evident from 
the development of a large variety of machines for calculating Fourier components. 
In the late nineteenth century and the first part of this century there were mechanical 
analysers working via a complicated system of gears and levers. Some of them can 
be seen in the Science Museum (London), including Henrici's Analyser of 1894, which 
is shown in the photograph. Descriptions of many analysers may be found in the 
literature (see Bibliography). The modern counterparts of these machines are the 
fast computer algorithms now available. 


* This function is quoted in E. C. Titchmarsh, The Theory of Functions 2nd ed., (Oxford University Press, 
1939), page 351. 
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Henrici's analyser Crown copyright, Science Muscum. 


The idea of representing a function f by an infinite series of known orthogonal 
functions ф, which form a basis of the function space was introduced in Unit M 201 20, 
Convergence and Bases. Such a series is called a generalized Fourier series. As was 
indicated in the Introduction to Unit M201 22, the most commonly used basis is the 
trigonometric basis which produces the (usual) Fourier series. We shall concentrate 
mainly on this basis, although we shall study first the more general case. 


6.2.1 Uniform and Mean Approximation 


The discussion of Section 6.1 suggests that in suitable circumstances the solution 
to a partial differential equation with appropriate boundary conditions may be 
expressed as an infinite series. The validity of a series solution depends on its con- 
vergence. We shall discuss three types of convergence in this section— pointwise, 
mean and uniform. The discussion of convergence will require some use of techniques 
from analysis. 


We shall meet the three types of convergence iri the next reading passage from W. 
However, we first describe uniform convergence in greater detail, since you have 
probably not met this concept before. 


The notion of uniform convergence is important in analysis and has played an 
important role in the development of the theory of generalized Fourier series, In 
this course, although it is not a main theme, a rough idea of it will be useful. 


SAQ 3 


Let s,(x) be the sum of the finite series 


xeR. 


By noting that this is a geometric series, show that 


1 


SX) = 1 + x? — —— 
S,(X) +x (tax 


(Solution on p. 27.) 
From this SAQ it follows that 


lim s(x) 21-x* хя 0. 


But s,(0) = 0, and so 


lim s,(0) ¥ lim lim s,(x) = 1. 


na x0 nso 


i.e. the two limiting processes as n ~ oo and x ~ 0 are not interchangeable, despite 
the fact that s,(x) is pointwise convergent for each xe R. 


This strange phenomenon was investigated in the nineteenth century by Stokes, 
Seidel and Weierstrass who showed that it cannot occur if the sequence of functions 
converges uniformly; in the above example the sequence s,(x), 52(x),... of partial 
sums does not converge uniformly in the neighbourhood of the origin, as we shall 
see soon. 


The important aspect of the uniform convergence of a sequence {f,} to f (which 
distinguishes it from mere pointwise convergence) is that given any £ > 0, 


f(x) — fied] < e for all n > №, and all x in the domain, 


where N, is independent of x. In the example above we showed that for х # 0 the 
sequence (s,(x)) converges pointwise to f(x) = 1 + x?. By the result of SAQ 3 


1 
|f = se) = Teer? 
so that given e > 0, 

If (x) — seo < £ 


only if n > N, where 


ing 


> =й ке еы ee сога 
Need In(1 + x?) 


This expression depends upon x, and is unbounded as x ~ 0. Thus the convergence 


of {s,} to f is not uniform. However, the series is uniformly convergent on the domain 
{x 2 хо) where хе > 0. 


We can represent the concept of uniform convergence diagrammatically as follows. 


ye Их] + к 


у= Их) 


= fx) E 


The equivalent situation for the sequence of functions s, defined in SAQ 3 which 
does not converge uniformly is shown in the next diagram. It is evident that the 
sequence is not uniformly convergent on any interval including the origin. 
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SAQ 4 
(i) Show that the sequence 
Alx) = п?х(1— xy 
has the pointwise limit, 


lim f(x) = 0, O<x<l. 

n ` 

(ii) Show that f,(x) has a local maximum at x = 1/(1 + n) which is unbounded as п 
becomes large. Hence show that the sequence { fp} does not converge uniformly. 


D 


(Solution on p. 28.) 


These examples show that uniform convergence is stronger than pointwise converg- 
ence, in the sense that the former implies the latter (this follows directly from the 
definition) but not vice versa. 


It is not immediately obvious from the above discussion why uniform convergence 
is important. There are three main reasons why it is, and we shall make use of them all 
in this unit. 


1 Continuity 


If a series of continuous functions is uniformly convergent in a given domain, then 
the sum is also a continuous function. 


Thus if u, (n = 1,2,...) are continuous functions with domain [a,b], and the series 


© 
Је = Xu 
п=1 
is uniformly convergent fora € z < b, then f is continuous on (a, b]. 


2 Integration 


If the series 


is uniformly convergent for a < z 
may be interchanged: 


[гаш = ИР eje - Y E 


qu] n=1 


S b then the order of integration and summation 


PDE 6.2.1 


That is, the series may be integrated term by term as though it were a finite series. 
(We shall see later that uniform convergence of the series is a sufficient but not necessary 
condition for this result to hold.) 


3 Differentiation 


If u(n = 1,2,...) are defined on the domain (а, b], and are such that uj(z) exist for 
all z in this domain, 


converges pointwise in the domain. and the series 


x 


У as) 


n=1 


is uniformly convergent for all z, then the function f given by 


© 


f= X u,(z) a&z&b, 


n=) D 


is differentiable and 


Thus, if the conditions of the theorem are satisfied the series may be differentiated 
term by term. 


That these results are important when expressing solutions of differential equations 
in terms of infinite series should not need emphasis. It should be borne in mind, 
however, that pointwise convergence is not sufficient for these results to be valid. 


Now we need some criteria indicating when a given series is uniformly convergent. 
There are two commonly used tests which we describe. We shall be using both of 
these tests later in this unit. 

CAUCHY'S TEST 


A necessary and sufficient condition for a sequence of functions и, with domain I 
to be uniformly convergent is that given any є > 0 there exists an integer N, inde- 
pendent of xe I such that 


lux) = u(x) < 5 


for all n, m > №, and all хє/. 


WEIERSTRASS'S M-TEST 


This is a sufficient though not necessary condition. If each term of the series 
го) = X ou) xel 


is positive and 
пх) < M,  WrelneZ*, 


where each M, is independent of x, and if 


а 


У, M, 


n=l 
is convergent, then the given series for f is uniformly convergent. 
We have presented these results in order to use them later. Proofs may be found in 
Kreider et al., Linear Analysis, Appendix 1-6. A more detailed discussion is given in 
Spivak, Calculus and Unit M231 15, Uniform Convergence. 
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READ W: Section 15. pages 70 to 72. 


Note 


W: page 71. line — 10 
More precisely. eigenfunctions have this property. 


General Comment 


In general. pointwise convergence does not imply convergence in the mean. This is 
demonstrated in the next SAQ. Conversely. as we have seen in Unit M207 20, mean 
convergence does not imply pointwise convergence; it follows that mean convergence 
does not imply uniform convergence. However, it may be shown that a uniformly 
convergent sequence of functions is convergent in the mean to the same limit. 


SAQ 5 
(i) Show that the sequence defined by 
1 
< — 
9 ш п+1 
с 1) -—— 
fox) = Шер nal Bu n 
1 
0 -«x&l 


has the pointwise limit 
lim f(x) 20 0Sx&l 


now 
Thus it is pointwise convergent in [0, 1]. 


(ii) Show also that 
1 
INCOLIS 
0 


for all n, so that the sequence does not converge in the mean to zero with respect 
to the weight function 1. 


Solution on p. 28.) 
SAQ 6 (Optional) 


Using the notation of И: page 72, line 8, suppose that f (x) ~ Y,,- .¢,@,(x) where the 
unctions ф, are orthonormal (i.e. orthogonal and such that 


b 
i $2p dx = 1). 


cote) = N к " 
Let sy(x) = ye 1C,0, (x). and put ty(x) = А. 15,@,(x), where the b, are arbitrary 
real numbers. Prove that : 


b е b 
a) L£(x) = зх) р(х) dx < || [/(х) — г(х]?р(х) dx 


(this verifies that the c,, as defined in W, provide the best-fitting approximation 
to fin the function space spanned by (di... Ф) 


(b) if N > M, 
b " b 
f І (F(x) — ssx) p(x) dx < | [/(х) — sy(x)P a(x) dx. 
(Solution on p. 29.) 
SAQ7 


W: page 72, Exercise 1 
(Solution on p. 29.) 
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+ 
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6.2.2 Parseval's Equation 


Having defined convergence, we can now investigate the conditions under which the 
Fourier series of a function converges. In the next reading passage some properties 
of the Fourier coefficients are obtained: the main results are conveniently stated as 
a theorem. 


THEOREM 


Let {ф,. ф›....) be orthonormal on [а,Ь] with respect to the weight function p. 
Suppose that # 


b 
| UC COT! p(x) dx 


exists, and that 


Са 


Јо) ~ X сф). i 


п=1 
Then the series 3; , c2 converges and satisfies 


b 


$ с р (р(х) dx. 


п=1 а 


This is Bessel’s inequality. Equality holds if and only if the Fourier series converges 
to fin the mean. The resulting equation is Parseval’s equation. 


The proof of a more general form of this theorem is contained in the first part of the 
next reading passage. 


READ W: Section 16, pages 73 to 75. 


Notes 

(i W: page 74, line 5 
We note here, without proof, that on the interval (— x, я) the set of functions 
{1, sin nx, cos nx:n = 1,2,...} forms a complete set. Note also that a complete 
set is just what was called a basis in Unit M 201 20. 


b 
A function f for which f f^ p dx exists is square integrable with weight function p. 
а 


(i) W: page 74, line 6 
The condition that f be continuous has been introduced to eliminate functions 
of the following type. Suppose that f(x) = 0 for all a < x < b except at one 
point where it is non-zero. Then 


b 
| /(х)ф„(Х)р(х) dx = 0 for all n. 


Thus this condition does not, on its own, imply that f = 0. The requirement 
that f be continuous allows the implication to follow through. 


(iii) W: page 75. line 6 | 
А series 


х 
У ay 


п=1 


of rea] numbers converges absolutely if the series 


È dad 


n=) 
converges. It follows that a convergent series each of whose terms is positive 
converges absolutely. It is proved in courses on analysis that if the terms of an 
absolutely convergent series are rearranged then the resulting series converges 
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to the same limit as the original series. and that two absolutely convergent series 
may be subtracted term by term. ‘ 
SAQS 
IV: page 76. Exercise 1 
(Solution on p. 29.) 
SAQ 9 
W: page 76, Exercise 2 


(Solution on p. 30.) 


SAQ 10 


For a string fixed at x = 0 and x = 1, and released from rest under the action of no 
body forces, the most general motion is of the form 


А 
y= У b,sin лях cos ппс, 


and the nth term of this series is called the nth harmonic of у. 


Using the expression for the energy in a uniform string of unit length derived in 
Unit 2, Classification and Characteristics (note (i) of Section 2.1.3), 


‘TP fay? ay\? 
РЕ Я i d ak | E ч 
Е = 3р Í 6 + с Е dx, 


where p is the line density of the string. show that Parseval's equation has the 
interpretation that the total energy of the motion is equal to the sum of the energies 
in each harmonic. 

You may assume that term-by-term differentiation of the series is valid. 


(Solution on p. 30.) 


6.2.3 Dini’s Test 


We now proceed to examine conditions for the pointwise convergence of trigono- 
metric Fourier series. The next reading passage shows that the Fourier series of the 
function / (with respect to the trigonometric basis {1, sin пх, cos ax ine Z*)) con- 
verges pointwise to f (x) provided 


Г If + т) - fool 
dt 
E а 


is finite (ie. exists). (This is Dini's test.) If the one-sided limits of the integrand do not 
exist as т ~ 0), then the integral has to be interpreted as 


lim Uc V + т) = f) Же Г If т) - fe) «| 


emo | Jag kal т 
Such an integral is called an improper integral. 


[The one-sided limits f(x + 0) and fix — 0) 
such that whenever т e (0, д) : 


Mix +) — f(x + OL «n 


exist, respectively, if Vn > 0, 36 > 0 


and 
[f(x т) - f(x — 0) < y. 


In Unit M201 22 these limits were denoted by f(x*) and f(x7) 
limits exist at x and are equal to Дх) then fis continuous at 
equal to f(x) then f has a jump discontinuity at x. 


. If both one-sided 
X; if they are not both 
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A piecewise continuous function on [a,b] is continuous at each point of [а,Ь] except 
for a finite number of jump discontinuities.] | 


SAQ 11 (Optional) 


Show that iff has a jump discontinuity at x such that f(x + 0) # fix = 0) then 
Г Ie т) fe 
dr 
-x kl 
does not exist. 


(Solution on p. 31.) 


In view of SAQ 11 it is desirable to be able to extend Dini's test to cover piecewise 
continuous functions, since these frequently occur in practice. In addition, such an 
extension is required at the end points +m. unless / (n) = f (—z). 


Dini's test is extended (in the next reading passage) to accommodate jump dis- 
continuities in f. At such points the Fourier series converges to 


if x + 0) + f(x — 0) . 
provided the modified Dini's test is satisfied. 


We shall also see that for suitable functions f, if f(x) exists for some xe[—, л] 
then Dini's test is satisfied at x. Extending this, we could derive the sufficiency 
condition (which you met in Section 22.2.2 of Unit M201 22) that the Fourier series 
of a piecewise smooth function converges pointwise to the average of the left and right 
limits. This result is the most important to remember; you should not spend too much 
time on the details in the next reading passage. 


READ W; Section 18, pages 77 to 80.. 


Notes 


(i) W: page 79, line 13 
We have not asked you to read W: Section 17 in which the Riemann- Lebesgue 
lemma is proved. The Riemann- Lebesgue lemma says that for suitable ortho- 
gonal sets of functions {фу} 


Г PnP 


lim туут 
t i Фә} 


provided 


b 
f [|р 


is finite (i.e. the integral exists). In our case 


= 0, 


a= –п, 
b=n 
р= 1, 
aye fix pu fi). 
sin jt 


dy(t) = sin (N + 1)т, 


and the set {sin(N + 4)t:N = 0,1.2,...} is orthogonal and "suitable". Since 


П 


i sin*(N + рил} = п! 


the condition that 


) 


fix * 0 - ft 


dt is finite 


sin ir 


is sufficient for the right-hand side of Equation (18.6) to approach zero. 


(i) W: page 79, Equation (18.7) 

The notation 
«x 

is used to denote that an expression is finite (or that a limit exists). 

(iii) W: page 79, line — 13 
The notion of Hólder continuity (more often called the Lipschitz condition 
of order g) is not important for this course. However, we remark that Hólder 
continuity implies ordinary continuity, but that the converse is not true (for 
example, the function x — (1 — x)7! is continuous on [0, 1) but is not Hélder 
continuous for апу о). We remark too that if x > 1, Hólder continuity implies 
differentiability. 


SAQ 12 


Why do the Fourier series of the functions f(x) = x", n being a positive integer, 
converge pointwise to f(x) for —z < х < л? What happens at x = +7? 


(Solution on p. 31.) 


SAQ 13 


Show that the Fourier series of the function defined on [— m, x] by 


Ш 
© 
| 
a 
IN 
E 
^ 
© 


fo) 


converges to f(x) everywhere except at x = 0, +z. What is the sum of the series at 
these points? 


(Solution on p. 32.) 


In W: Section 19 it is proved that if f is continuous on [—z, z], f (— n) = f (x) and f’ 
is square integrable on [ — n, x) (i.e. 


[^ 


is finite) then the trigonometric Fourier series of f converges uniformly to f. This result 
will be referred to in later reading passages: 


Weinberger then goes on to use this property of uniform convergence to prove that 
the orthogonal set of functions 


{1, sin nx, cos nxin = 1,2,...} 


2, 


is complete on the interval [—z z]. This is an important r А 
D^ esult which you 
remember. you should 


We have not asked you to read the proofs in W: Section 19 sin 
are rather more concerned in the applications of Fourier se 


Moreover, proofs applicable to a wider range of functions wil 
Sturm- Liouville Theory. 


ce, at this stage, we 
ries than the theory, 
l be given in Unit 13, 


A sufficiency condition for a Fourier series to conver 
represents is Dini's test. Its main applicabilit 
on [—2, x]. 


] 8€ pointwise to the function it 
y is to functions which are differentiable 


This concludes our brief analysis of the convergence of Fourier series, We have 1 
skated on the surface of the subject, a complete treatment of which would bie. | 


least опе full course and a lot of sophisticated analysis. is 
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6.3 GENERALIZED TRIGONOMETRIC SERIES 
6.3.0 Introduction 


We have seen in the previous section how to form trigonometric Fourier series using 
the complete orthogonal set 


{1 sin nx, cos nx:ine Z*). 


The importance of such series is evident when we consider the problems in Section 6.1 
which led to series expansions with respect to this basis. 


We shall now consider several generalizations which extend the use of trigonometric 
series. The first is the complex form which uses the basis {e""*:n e Z} for the vector 


space of functions [—x, x] — C. This space contains all the real-valued functions 
on [~2, л], and indeed the basis functions are closely related to the trigonometric 
functions by 


еї" = cos nx + isin nx nez. 
We shall also consider sine and cosine series, which make use of the fact that 
{sinnx:neZ*} and {cosnx:ne Zi) 


are complete orthogonal sets of functions on [0, л). Finally we shall investigate 
trigonometric series on an arbitrary domain [a, b). These last two topics constitute 
straight revision of material in Unit M201 22. 


6.3.1 Complex Fourier Series 


^ 
In this section we shall sketch some results on complex Fourier series, appealing 
when necessary to the relationship between the complex exponential function and 
the trigonometric functions given by Euler's formula 


e" = cos nx + isin nx 


(Unit M 100 29, Complex Numbers II). We obtain from this formula the results 


Ix E 
cos X = 5 (e* Te, 


€——À le — e7’) 
sinx 75i , 


e" -(-1" neZ. 


Our approach should be reminiscent of Unit M201 31, Fourier Transforms, where 
we used complex-valued functions to unify the handling of the sine and cosine 
transforms. In order to deal with complex-valued functions it is necessary to define 
a complex inner product. This is similar to a real inner product. Given a complex 
vector space V (ie. a vector space in which multiplication by complex scalars is 


allowed) and a mapping 
Vx ИС 


a complex inner product is defined on V, if for all a, b, ce V and 4, pe C the following 


axioms hold. 
1 а:а= 0=а = 0 


2 а:а>20 
3 а.Ь=Б-а 
4 а. (2 + uc) = Ха: + ware 


Note how the conjugate symmetry ахіот 3 ensures that а.а is real for all ae V. We 
are interested in the complex vector space whose elements are functions 


/{—т,л] — C: ' 
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As usual. we define two vectors a and b to be orthogonal if 
‘ 


a-b=0. 


SAQ 14 
Show that the functions 
xe) neZ 


are orthogonal with respect to the inner product 
f-g= | [(х)в(х) dx. 
E 


(You need not verify that the formula defines a complex inner product.) 
(Solution on p. 32.) 


The set of orthogonal functions {e'*} is equivalent to the trigonometric basis because 
we have chosen a complex inner product which reduces to our usual real inner product 
when fand g are real-valued functions, and Euler's formula represents e*'** as linear 
combinations of cos nx and sin nx, and vice-versa. Consequently, it may be shown 
that Ге" :n& Z} is complete. We could follow through our discussion of Fourier 
series for this set of functions, and we would obtain 


E 


Ло) X ce". 


п=-ж 


$АО 15 
Show that 
Ir 


Cy = = 
2n Ju» 


e^ Pf (x) dx. 


n 


(Solution on p. 32.) 


The advantages of the complex Fourier series are many. Firstly, it is less space 
consuming to write down the complex form instead of the real form, and the two 
calculations for the a, and b, are combined in one calculation for the С, 50 that only 
half of the work is required; we shall see an example of this in the solution to SAQ 18. 
Secondly, under the operations of differentiation and integration e" is easier to 
handle than either sin nx or cos nx. Thus, 


qP 
IL nx) = (—lyn?sinnx if p= 24 
x 
= (-l'mcosunx if р=20+1 
а? z А 
Jyp (COS nx) = (—1)°созпх if p= 2q¢ 


= (= 1) sinnx if p=2q-1 


È oes cob dd 
dx? 

with similar results for integration. The usefulness of this simplification where 
dealing with differential equations need hardly be mentioned. 


There are other reasons why the complex form of F 
theory of real analysis, of which Fourier series is 
viewed from the more general complex variable t 
tial function is a fundamental function. 


ourier series is used. Most of the 
a part, is easier to understand when 
heory and in this theory the exponen- 


In obtaining these simplifications we have introd 
numbers; but when one is used to their manip 
countered when handling sines and cosines, whi 


their 
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uced the inconvenience of complex 
ulations they avoid the pitfalls en- 


sspe mudileda ch are prone, at the least, to getting 


540 16 


If 


and 


- do z К 
f(x) ~ >t Yt, COS nx + Б, sin nx) 
= 1 


show that 


Cy = Ha, — iby) 50 
п & 
Coa = Ha, + ib,) 


ё 1 
Co = 209. 


(Solution on p. 32.) 


540 
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Show that if f is a real-valued function then 


Oy = Cun. 


(Solution on p. 33.) 


6.3.2 Sine and Cosine Series 


The 
even 


gist of this section is that an even (odd) function can be expressed as a sum of 
(odd) functions. Since cos is an even function and sin is an odd function the 


Fourier series for an even (odd) function contains only cosines (sines). 


REA 


Note 


SAQ 


D W: Section 20, pages 87 to 88. 


W: page 87. lines —6 and —5 

The reference to Parseval's equation is unnecessary. Since |l. sin пх, cos nx in = 
1,2....} is complete (as proved іп W^ Section 19) any suitably behaved function 
can be expressed in terms of these functions with Fourier components given by 
Equations (18.1) in W: Page 78. If the function is odd, then а, = 0 for all n. and 
it follows that any suitably behaved odd function may be expressed in terms of 
sin nx (1 = 1.2...) only. Thus fsin nx :n = 1. 2....j is complete on the space 
of odd functions with domain [—2, x]. Since any function with domain (0. я] 
can be extended as an odd function on [ — x. л]. (sin nx} is complete on the space 


of (square integrable) functions with domain [0, л]. 


18 


W: page 88. Exercise 1 


(Solution on p. 33.) 
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6.3.3 Change of Scale 


READ W: Section 21, pages 88 to 92. 


SAQ 19 

W: page 92, Exercise 3 
(Solution on p. 34.) 
SAQ 20 

№: page 92, Exercise 5 


(Solution on p. 35.) 
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6.4 APPLICATIONS 
6.4.0 Introduction 


We now proceed to apply the methods discussed in the previous sections to problems 
involving partial differential equations. Many of these applications were introduced 
in Unit M201 32. where both the heat equation and Laplace's equation were solved 
fora variety of boundary conditions. The treatment presented here is slightly different, 
much more effort being devoted to rigour. 


The problems which concern us involve a partial diflerential equation in some 
domain D with a continuous boundary condition on its boundary C. In the case of 
the heat equation, С may be taken as the physical boundary for г > 0 together with 
the whole initial line. A solution of the problem must satisfy the equation and 
"boundary" conditions and be continuous on DUC. 


In W: Sections 22 and 23 the following method is used. Firstly, the method of 
separation of variables yields a Fourier series solution which formally satisfies the 
differential equation and subsidiary conditions: it is still necessary to prove that this 
series does indeed converge to a solution. Next. Weierst M-test (Section 6.2.1) 
is used to show that the series is uniformly convergent on part of D u C. This ensures 
that the sum u of the series is continuous in the relevant part of D ù C. Next, itis shown 
that the series obtained by differentiating term by term at any point in D is uniformly 
convergent in some subset of D. so that u really satisfies the differential equation. 
Finally, Cauchy's test (Section 6.2.1) and the maximum principle are used to ensure 
uniform convergence throughout DUC (provided some extra conditions are 
satisfied) so that и is continuous everywhere. 


Note that in Unit 3 we have already proved, using the maximum principle, that 
suitable subsidiary conditions for the heat equation and Laplace's equation yield 
uniqueness of the solution and its continuity with respect to the data. This section proves 
the existence of solutions to selected problems of this type and thus shows that such 
problems are properly posed. This is as we would expect, since ihe problems are based 
on realistic physical situations. 


6.41 The Heat Equation 


Our first example is the one-dimensional heat equation with homogeneous boundary 
conditions. 


READ W: Section 22, pages 92 to 95. 


Notes 
(i) W: page 93, lines 15 to 17 


The statement that (b,] is uniformly bounded means that there exists сє R, inde- 
pendent of n, such that, for all п, |р, < c. 


The uniform bound is obtained by using the result 


b " 
f g(x) dx 


which is the integral equivalent of the triangle rule for finite sums. 


b 
X | Igi dx 


а 


[п uso ab S ded + deal + lanl- 


Then since |sin nx| < 1. | fx) sin пх| S |/(Х)| and 


2 
< | LEGO dx. 
T Jo 


D, 


9 рт 
Б, = ef fix) sin nx dx 
Tdo 
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{н} We page 93. lines —11 to -4 "——" 
You need not verify that these series converge for г > 0. However. this sai ү 
is obvious that the convergence is uniform in x. (To show that it is uniform for 

H ; Ж үза amt. - the 
LÈ to. we require the M-test with M, = exp( — n^Kto) and similarly for the 
other series.) 

(ili) H^: page 94, lines 10 to 12 : . 

This follows because the Fourier sine series of fis uniformly convergent. 


SAQ 21 


Find the analytical solution of 


ыы c O0cxc«Lt-90 
Ct сх” 
u(0, 1) = (1.0) = 0 г> 0, 
2х 0<х<} 
шыш e =x) i«x«L 


We discussed the numerical solution to this problem in Unit 5, Initial Value Problems 
(S: pages 11 to 16). 


(Solution on p. 35.) 


SAQ 22 


Write down the heat equation and boundary conditions for a thin conducting ring 
or armlet of circumference 2/ if the initial temperature is sin?(nx/21), where x is the 
distance around the circumference from some fixed point. Determine the subsequent 
temperature distribution. 


This problem is especially interesting as it was one of the first to which Fourier 
applied his mathematical theory of heat, and for which the results of his mathematical 
investigations were compared with experiment. For the original theory see Fourier, 
The Analytical Theory of Heat. 


(Solution on p. 36.) 


6.4.2 Laplace's Equation 


READ W: Section 23, pages 95 to 98. 

шы r + = шы IL == ae. — —- TE 
If you are short of time you may omit the discussion of error bounds, W: page 97, 
line — 11 to page 98, line 7. 


Note 


W page 98, final paragraph 

Section 25 is not included in this course. 
As an illustration of this method of solving Laplace's equation we consider the 
problem discussed in Unit 3. that is, the calculation of the total rate of steady flow, 


Q. ofa viscous fluid through a pipe of square cross section. In Section 3.2.3 we showed ` 
that 


Q 
$55 <= <0, 
0555 < E, < 0.667, 


where 2a is the length of a side of the square and К a constant (defined in Section 3.2.1). 
Here we determine the exact value of Q which is given by 


Q -Í f wdx dy 


кә 
е) 
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Where w is the velocity of the fluid and satisfies the equation 
Viw = —К, 


with w = 0 on the boundary. This equation is not quite of the form we require, but 
on defining 


dx. у) = wt ee +y? — 2a?) 
we find that 
Vie =0 
and that 
ф(х,у) = р + у? — 2a?) on the boundary. 


The choice of ¢ satisfying V?¢ = 0 is not unique; any function of the form 
к. › 
(x, y) эз + ris + у) + ax + Ву уху д 


where g, fj. y and à are arbitrary constants will do (since the second term goes to k 
2 vis 

under the operator V^, and the remaining terms go to zero). The reason for our 

choice will soon be apparent. 


Although now the equation to be solved is simpler the boundary conditions are more 
complicated. In order to solve this problem we use the technique of W: page 98 
(also discussed in W: page 64) and put ф equal to ф, + ф› + фз + ġa. each of which 
satisfies V7; = 0 and is nonzero on only one side of the square. This apparently 
necessitates solving four equations. But the symmetry of the problem shows that only 
$, need be found; all the others can be obtained by either changing signs, or inter- 
changing x and y. Thus we need to solve 


Vie, =0 X,ye(—a,a), 
koz 2 

ф(х, a) = gh = а?) xe[-a,a], 

$,(x,—a)20 xe[-a,a], 


pilay) = ф(—ау)=0 rel-ad. 


Now it is clear how the arbitrary constants of ф were chosen. With our choice of 
constants, ф(х, y) = 0 at the corners of the square so that each ф, is continuous on 
the boundary. With any other choice of constants $; would not be continuous on 
the boundary, and discontinuities are best avoided as the Fourier series of dis- 
continuous functions are not uniformly convergent over the, whole interval, and 
also converge more slowly, creating a new computational inconvenience. 


The solution to this problem is now straightforward but lengthy. By separating the 
variables and using the boundary conditions along x = ta and у = —a we find that 


= . nn(y +a). mix + a) 
х,у) = „sinh - — sin . 
Фк) = Y e, 2a 2a 


п=1 Es 


Aty =a, 


E x . . nn(x + a) 
(x? — a?) = У c, sinh um SIDE 


n=1 ^ 


xe[-a.a]. 


Q(x. а) = 


REX 


so that 


1° k . nn(x + a) 
c,sinh nz = — - (x? — а?) ѕіп ср Cl 
ај _.4 2а 


dark 
= CY 1). 


mmu 
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Since the boundary condition is continuous. the serics for $, is uniformly convergent 
in [а.а] x [—a.a} and the differentiated series are uniformly convergent in 
(—u.d) X (аа). as in the reading passage. 


The total flux is now given by 


Q= Ko l^ Gr) dux y) "— d dx dy 


4kat 
=4 ii f pixy) dx dy + EB 


р 


+4 St sin —— 


3 n=1 -а 


ке]; _ 256 cosh on +1)л—1 | 


п^ „со (2m + 1) ѕіпһ (2m + Iz 


gka? ©... "е + а) dx EL п y+ а) 
2а 


Using the first three terms, this gives 
О = 0.562ka* 


to three places of decimals. 


SAQ 23 
W: page 99. Exercise 9 
(Solution on p. 37.) 


Laplace's equation may be solved for a circular boundary by methods similar to 
those of this section. The discussion appears in W: Section 24 which we shall cover 
in Unit 10, Green's Functions 11. 


PDE 6.5 
6.5 SUMMARY 


In this unit we have considered two main topics, both of which were first introduced 
in М201. 


We showed how some partial differential equations could be solved by separation of 
variables, We saw that equations which permitted this treatment led to Fourier series; 
and we found that certain boundary value problems and initial-boundary value 
problems yielded solutions under this treatment. 


In considering the theory of Fourier series we introduced the notions of uniform 
convergence, convergence in the mean and pointwise convergence: the ideas of ortho- 
gonal, orthonormal and complete sets of functions were also used. A new result, 
Dinis Test, was obtained as a sufficient condition for the pointwise convergence of a 
trigonometric series. We also obtained a general form of Parseval’s equation which 
implies that a Fourier series (with respect to any complete orthogonal set of functions) 
can be integrated term by term to yield a pointwise convergent series. 


These techniques were applied to a few physical examples, usingsome general theorems 
on uniform convergence and a result (which we did not prove) concerning sufficient 
conditions for a trigonometric Fourier series to converge uniformly. 


A brief discussion illustrated the use of complex exponential functions as a complete 
set of orthogonal functions which yield the complex generalization of trigonometric 
Fourier series. 


m 
D 
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6.6 FURTHER SELF-ASSESSMENT QUESTIONS 


SAQ 24 


Show that, forz > =}. 


lim i x^ sin nx dx = 0. 
0 


nc 


(Solution on p. 38.) 


SAQ 25 
The function fis odd and for 0 < x < л, f(x) = cos x. Obtain the Fourier series for 


fin the form 


А 8 5 msi X 
fiie p am = 1 


т= 1 


(Solution on p. 38.) 


SAQ 26 


A function f with period 2л is equal to x? for — z € x € л. Show that 


x п? be 
fw=>+4 У 


n=l 


COS AT COS HX 
Eee ees xeR. 
t^ 


Deduce that 


x] 2 


(Solution on p. 39.) 


SAQ 27 
Is the function 
fix > |х| xe(—7,0) o (0,2) 


piecewise smooth for ж > 0? Show that if 0 < x < 1 then the Fourier series for f 
converges pointwise to f for x # 0. ` 


(Solution on p. 39.) 


26 
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6.7 SOLUTIONS TO SELF-ASSESSMENT QUESTIONS 


Solution to SAQ I 

Substituting y = X (x)T(r) into the equation we obtain 
PE Gu d іХ 
Xt) (0 = TO ^ 1X о) 

which becomes, on rearrangement, 


1 ФТ 1 d[ dX 
cus ) mex 
T(t) d? X(x)dx \ dx 


(x) 


Since the left-hand side depends on г only, and the right-hand side is independent 
of t, each side must be equal to a constant, а say. Thus the required equations are 


dT T=0 
dt? TUN 
dX dX 
Ф? tdg 7-9 


Solution to SAQ 2 


Since the solution possesses radial symmetry the most suitable coordinates are polar 
coordinates (r, 0) in which the solution must be independent of the polar angle 0. 
The equation is thus 


И 
Since ф has axial symmetry, 2ф/20 = 0 and we may write 


Pb dọ _ : * 
(rtt Ocr«l, (*) 


where $ depends on r only. 


We look for a solution of the form ф = :* where « is some constant; on substituting 
into the equation we find that 


(02 — 1y* = 0. 
Since i? 4 0, a = +1, and the solution space contains r^! and ғ. Since the equation 
(*) is normal} and second-order, its solution space is two-dimensional. So the general 
solution is 

ф = Afr + Br. 
But the solution is bounded in (0, 1), and this can only be so if A = 0: also, at r = 1 
we have ф = B = 1, giving 

par 


as the required solution. 


Solution to SAQ 3 
The sum of a geometric series is given by 
onl 
ay =a 
pro 
In this case a = x?, r = 1/(1 + x^); thus 
(x)= 14x? 
SAN) х= 1 x ае 
с (1+ х2)": 
+A differential equation is normal on the domain / if its leading coefficient is nonzero throughout / 
(Unit M201 9, Differential Equations 11) 


Р ОЕ SAQ 


Solution to SAQ 4 


‘ 
ü) For x = 0.1 we have /(0) = 10) = 0 for cach ne Z*.so that 


lim f(x) = 0 х= 0.1. 


Let£z 1 — x. ПО < x « 1 then0 < č < I. so that 
0 « f(x) = хп?" « n^i" O<x<1 
and, since |2] < 1. 


lims?" = 0. 
PE 
(To verify this statement. note that 0 < ët < 1 so that we may choose an integer 
N such that ` 
E N 
S e Es 
N+1 


> 


For п > N the sequence 32" is decreasing. i.e. 


(п + PEt) < mP < NEY = k, say. 
Thus 
] Р ‚К 
0 € lim wë" € lim - =0.) 
nx nwa П 
Hence 


lim f(x) = 0 O<x <i. 


(i) The stationary points of f(x) are given by the roots of f(x) = 0, i.e. 
Ax) = aU х) хн) — xy7! = 0. 


Thus f(x) has a stationary point at x = 1/(1 + n) which can be seen to be a 
maximum by inspection. At this point 


Е п? 1 

A0) 7 гаг A 
Since 

е li 

lim |1 4 -| =e 

ane n 


(Unit M100 7, Sequences and Limits 1), УД + п) behaves as n and so is 
unbounded. 
Solution to SAQ 5 
(i) For pointwise convergence to 0 we require, for each x in [0, 1], that: 
given any c > 0.3N such that LX < efor all n» N. 


Now. given any x in (0, 1] we can always find an N » 1 such that x > 1/N 
(take N to be the first integer > 1/x). Then, for n > N f(x) = 0 forO<x <1, 
It follows that 


lim f(x) = 0 O<x <1, 


апа from the definition f,(0) = 0 Vn, so that the result follows. 
(ii) We have 


1 im 
| LG? dx = | a(n + 1)4х 
о lent 
= n(n + 1) [ m ттт =] 
n nl 
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This SAQ illustr. 


ates the point that pointwis r ; 
; Wise convergence does i BOR VORNEHAR: 
in the mean. gi oes not imply convergence 


Solution to SAQ 6 


(a) Expanding the right-hand side and completing the square as in W: page 71, and 


h 
remembering that [| Фар dx is unity since the Pa are orthonormal, we obtain 
a 


" 
| LPO) = rep! p(x) dx 


х b 2 b х 
XD - | ГФ,р as} + | fads = X { 
" 


ph 


) 


N b f, 
z (b, = e | f2pdx — У c by the definition of c,. 


п=1 


ll 


net 


Гфьр ax} 


1 


Thus, 


RN b х, 
| UG) = ss] р(х) dx = | J pdx- Y e, 


п=1 


ѕо that 


b h N 
[ Uf) = sso) plx) dx = | Lf) = орх) ах Y, (b, — cy. 


n=l 


and the result follows. 
(b) This result follows directly from part (a) by putting 
b, = с 
Б, = 0 ї= М+ 1... №, 


ї = 1,2,...М, 


so that гух) = у(х). 


Solution го SAQ 7 
We need to show that 
л 
i sin nx sin mx dx = 0 nám nmez'.. 
0 


The integral may be written in the form 


1" 
з] (соз(л — m)x — cos(n + m)x) dx 
o 


nám 


2 п = т п+т fo 


1 Е = т)х  sin(n + ды, 


= 0. 


For л = т the integral is 


л 
‚з 
| sin? nx dx = 
0 


Solution to SAQ 8 


Nia 


The Fourier sine series for 


Д(х) = 1 О<х<я 


Јо) ~ У d,sin ax, 


п=1 


39 
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where the Fourier coeficients are given by 


ex 1 n 
pm nxdx -| —cosnx 
n o 


i sin? nx dx п/2 
о 
| 0 neven 
4 
B — n odd. 
zn 


Thus 


e [0. zr]. 


If the set of functions {sin nx} is complete on (0, л] we may use the result of Equation 
(16.8) in W: page 75 to integrate term by term. obtaining 


* Ix’ if sin(2k — sin(2k — Dx Ix 
v х 
fee eS Ram 0 2k — 1 т-р" 


tale cos(2k — 1)х . 
x= x2. Oke x € [0. л]. 


or 


The equality here denotes pointwise convergence. 
Solution to SAQ 9 


Parseval's equation gives 


т 16 5 1 z 
Ix = = in?(2k хах, 
f dx 5 У Gk Fh sin*(2k — 1)x dx 


о п? үс 


Ия a, Ok EET 


Substitution of л into the series for x in the solution to SAQ 8 yields the same result: 


Ё, Rk- 8 


Solution to SAQ 10 


The energy of a uniform string of unit length is given as 


эте к 


The motion is of the form 


A 


= Y b, sin nax cos плет OSx< 1120, 


n=l 
and the nth harmonic is 
„ = b,sin mx cos naci. 


The energy Е, in this harmonic is obtained simply by substituting y 


in the above 
expression for the energy: " | = 


30 
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md 2 2 Tm А, 
E, = уры (ппс) || (sin? nax sin? nnet + cos? nax cos? net] dx 
о 


1 1 

are ‹ 2 " ; 

= yp(nnch,)? since [ sin? nax dx =| cos! nax dx = 1. 
0 [uU 


Now, term-by-term differentiation yields the Fourier series 


— У ancb, sin nax sin naet 
п= 1 


апа 


-= Y nab, cos ллх cos пле. 
n=] 


and applying Parseval’s equation, we obtain 


(i 


ж 1 
i : 
У (nach,)? sin? imet | sin? nax dx 
п=1 ° 

and 


2 ЕА 1 
dx = У (nnb,)* cos? пле | cos? плх dx. 
n=l о 


The total energy is then 


E 


I 


H 1 

1 a е 

4р D (nzcb,) since | sin? nzx dx = | cos? nax dx = $ 
n=l 0 o 


Il 


Y E. 


п= 1 
Thus the total energy of a vibrating string (or more generally a vibrating system) is 
the sum ofthe energies in each harmonic. 


Solution to SAQ 11 

Let = if (x + 0) — f(x — O); suppose that 
|f (x + 0) — /(х)| > 2. 

(If not, then 
х — 0) — f(x) > 29. 

and a similar argument will follow.) Now 20 > 0 such that 
Vx e (0. д) [fs +t) fix +0) < 1. 

Hence, 


Vre(0,2) Axt — fi)» n > 0. 


Therefore 
Siriy — fix $1 | 
lim Ше e n - fe It > y lim | -dr = 4 lim In(à/z) 
1007 Jp |т| є=0* de T 2-07 
which does not exist. Hence. 
Г еж 0) fe 
does not exist. F 


Solution to SAQ 12 


Because the functions f(x) = X" are differentiable for —z < x « x and absolutely 


integrable. 

For n even f(x) = f(—2). with the consequence that the Fourier series converges 
pointwise 10 fix)for -x €x € m. If n is odd f(z) # f/(—70. and the Fourier series 
converges to zero for x = їл. since }[ f(a) + (77) = 0. 
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Solution to SAQ 13 


i 
The function is piecewise smooth and so at each interior point converges to 


x40 E 0 -n«x«0 
ILU ы: а 1 О<х<л 
я lo x20 


whilst at x = +7 it converges 10 


? 2? 


Solution to SAQ 14 


Stan + 0) + fir- 0) _ 1 


The solution to this SAQ is similar to that of SAQ 7. We need to evaluate the integral 


[= | qi gs n.meZ. 
-7 


For n = m clearly I = 27. since e? = 1. Forn # m 


— mm 
i(n — n) R 


1 


РЕ (очи) отіп ту 
i(n — m) 
„ілби ml 

= a (1 — етт) = 0), 
i(n—nm) 


since e??? = | when p is an integer. 


Solution to SAQ 15 


The main difference between our case and the discussion in W : Section 15 is that the 
exponential functions are orthogonal with respect to the complex inner product 
defined in SAQ 14. (The complex inner product is more general than the real one 
defined in H^) Following through the analysis of W: Section 15 with our complex 
inner product then gives 


Í Sp dx 


c = SEENE 


NOE 


In our case ф(х) = e", so that 


| | Гох)е dx 


1 i i 
ye ————— = x Гое dx, 
| dx n: 


Solution to SAQ 16 


We have 


Cy = 5 | JF Ge" dx nez (SAQ 15) 
and 
Lf* , 
a, = aT f (x) cos nx dx nez; 


(Equation (18.1) in РУ: page 78). 


Ш 


LI y 3 
м = | Д\х) sin nx dx nez* 
2-а 
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Thus. for 1 > 0 


Ya, — ib) = 


f Г(х)}(соѕ nx — isin nx) dx 


Lo В 
= = Г Seve Dix dx 


= [^ 
and 
Ma, + ib) = LT p C+ isi 
200, п) 73x ul (x)(cos nx + isin nx) dx 
i T* 2204 
= 2x f J (х)е" dx 
= Cn: 
Also 


im 
м = y | rends 


Co- 
Solution to SAQ 17 


Since 


Тир : 
e =з- Гете dx, 


we have + 


Te dx 


л 
-r 


Í 


vena 


Гохе” ах since (x) is real 


Solution to SAQ 18 

Both of the required series may be found in one calculation if the function 
Д(х)= х О<х<& п 

is extended into 
[77.2] 


in the correct manner. For the sine series we note that sin nx is odd and so we extend 
[so that it is odd, i.e. 


Дх) = х —л<х<&л. 
The cosine series is even and so we extend f'so that it is even. i.e. 
Jax) = |х| —-n &xXm. 


If we now find the Fourier series of f(x) and f(x) on —z € x < л, they will be the 
required sine and cosine series for 


Го) =x О<х& л. 


ре 
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Thus we have 


X ^ Y b,sinnx 


n-1 
where 
2p 
Б = | x sin ax dx, 
To 
and 
o s 
lx] oe, а, cos nx 
ж nz] 
where 
X cos nx dx. 
Writing 


Ca = d, + ib, = al хе!" dx, 
and integrating by parts, we find that 


e= Je esc »| n> 0, 


n 


in 
and so 
do = п, 
Thus 
sei 
хе -2 У ——sinnx on [0л] 
ai 1 
and 
n 44 x 
“~i Taa к-ту Оп ad 


Solution to SAQ 19 


Let x and г be the space and time coordinates for the first bar. The heat equation 
describing its temperature u(x, г) is 


where i(X, 1) = u(x,t). This is just the equation satisficd b 
second bar with space and time coordinates denoted by 
therefore take as 


y the temperature of the 
X and 7, which we may 


au(x, t) + b 


34 


PDE 6 SAQ 19 21 


(using the hint), where a and b are arbitrary constants. To determine these we notc 
that att = 7 = 0 we have 


ullo) = Tp, 
ап(31.0) + b = Ty. 
and at x = X = 0 we have 
u(0. 1) = T, 1 
а1(0. 1) + b = T, І 
Since R(X, 1) = u(x. t). we obtain 
Ty = aTy +b and T, — aT, +b, 


giving 


The temperature of the centre of the second bar at time г is given by 
аши) + b = af PR) + b 
(To = TUPA + (ToT, — ToT), 
Т - T, 


Solution to SAQ 20 


Consider the transformation x = gX, t = fft, under which the equation becomes 


where ü(X,T) = u(xX, BT). Thus on choosing f! = 1/а and c? f*/g? = 1. or x = c/a 
we obtain the required result. 

Solution to SAQ 21 

We require the solution of 


A a? 
си он 
= 5 = 0 Ocx«l г> 0, 
ð ex 


with the conditions 


(0,0) = u(1,7) = 0 t>0, 
2x o<x<} 
s0-6 5 ES x] 


We consider a solution of the form u(x,t) = X(x)T(t). which on substitution into 
the heat equation gives 
1dT 1 d?X 
TETY 
where the constant —g? has been taken to be negative for reasons which will be clear 
shortly. 


2 
—-——W, 


The solutions are then 
T(t) = Ает, 
X(x) = Bsinax + C cos 2х. 
The boundary conditions at x = 0. 1 give: 
X(02C-0 
X(1) = Bsing =0 
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so that x = na(n = 1.2... ) and the general solution is 


щх = У A, sin imxe 


ned 


(If we had not chosen —z? as the constant, we would have obtained the trivial 
solution и = 0.) In particular when т = 0 


щх.0) = У 4, sin nax. 


aed 


But for 0 € x € 1. мх. 0) is given, and we have 


1 
A, = af u(x, O) sin aax dx 
0 


i 1 
4 x sin nzx dx + af (1 — x)sin nax dx. 
0° Н 


Substituting 1 — x for x in the second integral we find that 
A, = 4(1 — i— 10) | х sin nax dx. 
0 
since cos пж = (— 1)". It follows that if is even A, = 0. 


The integral can be integrated by parts to give 


A 8sin((2n +1 
Ber A 


and so 


S. 1 ‚ (Оп + ж. 22 
sin! LS т sin(2n + I)nx exp( — (2n + 1)?°л?г), 


which has the required convergence properties since x +> u(x, 0) is continuous, 


1 1 
u(0,0) = u(1.0) = 0, and | 4dx + | 44х is finite. 


о 
Solution to SAQ 22 


If x measures the distance around the circumference, the heat equation is approximately 
(for a thin ring) : 


where the conductivity has been put equal to unity, and u(x, г) gives the temperature. 


When x = 2/ we return to the point x = 0. so w(x, 1) must be periodic in x with period 
3l: 


u(x.t) = u(x + 2L) t20. . 

As in the solution to SAQ 21 we find that a separated solution is given by 
e^" А їп ах + В cos ux) 

and since this must be periodic in x, it is necessary that 


"n 
@=-— 


l 


Thus the general solution is 


ux.) = Ag + Y A, cos T XxX + В, sin exp| - ЕЙ 


п=1 
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and at time t = 0 this is 


u(x.0) = Ag + Y (a. cos = + В, sin = 


n= 


1 


sa TE. 
sm — 
5 


1 
io 
| 
с 
© 
a 
EI 
4 
st 


By equating coefficients we obtain 


1 
u(x, t) = 5 |: — eT Dt cos PX 


Solution to SAQ 23 


Using the method of separation of variables and writing u(x. y) = X(x)Y(y) the 
equation becomes 


гах " 14Х 1Y 
Хайй туы Үн С 


with 
Y(0) = Ү(п) = 0 
and 
X(0) = 0. 
As usual we write 
X" + X'-14X = 0 
and 
YAY = 0. 
The boundary conditions on Y yield the eigenvalues 
ice neZ*. 


and 
Y) = sin лу. 


Let X, be the general solution of 
Ху + X,- PX, = 0. 
Then we set 
з 
u(x.y) = У X,(x)sin ny: 
п= 1 
the boundary condition 
u(r. y) = sin y 
gives, formally. 
sin v E: X (a) sin ny. 


n=l 
where the Fourier coefficients come out as 
Xim) = 1. Хт) = 0 nl. 
Solving the boundary value problems for X, we obtain 
day sinh (/5x/2) 
Хх eset oe E p 
sinh (4, 57/2) 


0 n> 0. 


Ш 
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Thus, 
9 /§x/2 
u(x, y) = etm Minh Co sin y 
sinh (/ 57/2) 


The validity of the solution is guaranteed. since it is clearly continuous on [0, z] x [0, я] 
and sufficiently differentiable in (0, x) x (0, x). 
Solution to SAQ 24 


In Equation (16.4), W: page 73, we set p = 1, ф(х) = sin nx. a = 0, b = п: thus we 
we have, since [sin ux :né Z* | is complete on [0, л], 


[f J (x) sin nx i g 
o 


lim -—— ~ =0 


sin? nx dx 
o 


it | J (XY! dx is finite. Putting f(x) = x* we find that 


0 
m ^ x 1 z 
i f(x) x= zs i 


which is finite ifa > — 1, and since 


E boy л 
sin хх = 
о 2 


it is evident that 


л. 


lim | x*sinnxdx = 0 fora > —i. 


пто Jo 


Solution to SAQ 25 


Since f is odd its Fourier series may be expressed in the form 


© 


Го) ~ Y b,sin nx, 
=1 


n= 


where 
"EL 
b = zf sin nx cos x dx. 


This may be written in the form 


AR. 
b, HI (sin(n + 1)x + sin(n — 1)x) dx 
0 


1 cos(n + 1)x i cos(n — 1)х | " | 
п nl п 1 o us 


1 1 1 
aait hu ey quor mes Nt 
VD 1) [ * | 
Thus b, = O if n is odd (41) and ifn = 2m 

8m 


b, = — 
^ — m(Am? — 1) 
Also, 
Эре Wont edt 
b, -i[ sin xcosx dx = = zm =0. 
Tdo л 2 jo 
Thus 


8 = msin2mx 
(4m? — 1) 
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Solution to SAQ 26 


Since fis even its Fourier series will contain only cosines, and the coefficients will be 
5 


л 

z 2 

a, == X* cos nx dx n=1,2,... 
T Jo 


2 а 
== ref х?е" dx, 


ш o 


where Re means the real part of. The integral is best integrated by parts: 


я 2xe"" дот] 
=. 


| xwe dy = 
о 


so that 


а = 


Also. 
2TT a 2л? 
а =- | х°4х = d 


The Fourier series converges to /'(x) at each point, so that 


л? = 
4 COS ПТ COS nx 
Јо) = 4ў ——3 


nel n 
Since f(x) = f (— n) the Fourier series converges to л? at x = л, and we have 


cos? пт 


Solution to SAQ 27 
The derived function of f is 

f'xea(-xy! xe(—m, 0) 

= 7! : x € (0, л). 

As x ~0* or 07, f' is unbounded since x > 0 and so the limits 

70 + 0) and /'(0— 0) 
do not exist, i.e. fis not piecewise smooth. 
Nevertheless, we have seen in W: page 79 that it is sufficient if Dini's test is satisfied; 
this will be the case at each point where f is differentiable. provided f is absolutely 


integrable. Now 


pn 


l.. 


IAD dt = | Ju? dt 


т 
=2 lim [| t77 dt since |f] is an even function 
e~o’ 


1+1) 
= 2 lim | —-- 
s~o 11-00, 


is 
=] Cg! provided 0 « a < I. 
== 
Thus, if 0 < x < | then the Fourier series of f converges 1o f at each point in 
(72.0) o (0. х). 
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Set Books 


G. D. Smith, Numerical Solution of Partial Differentiation Equations (Oxford, 1971). 
Н. F. Weinberger, A First Course in Partial Differential Equations (Blaisdell, 1965). 


It is essential to have these books: the course is based on them and will not make 
sense without them. They are referred to in the text as S and W respectively. 


Unit 7 is not based on either set book. 


Conventions 


Before working through this text make sure you have read A Guide to the Course: 
Partial Differential Equations of Applied Mathematics. References to Open University 
courses in mathematics take the form: 


Unit M100 13. Integration П for the Mathematics Foundation Course. 


Unir M201 23, The Ware Equation for the Linear Mathematics Course. 
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7.0 INTRODUCTION 


This unit is essentially a case study 


using some of the techniques developed earli 
. - . B t а e 
in the course. In it we investig. i п 


| 1 | ate the modelling of the behaviour of an overhead wire 
when a high-speed electric train passes. We start with a brief history of the develop- 
ment of railways and their electrification. 


7.0.1 History of the Electric Train 


The earliest railways in Britain are dated at about 1550; by the end of the sixteenth 
century horse-drawn trains were common, particularly in the North-East coal 
mining area, and by the eighteenth century horse-drawn trains were found all over 
Wales, Scotland and England. The first steam locomotive. built by Trevithick, ran 
in 1804 but was too heavy for the track which kept breaking; and in 1823 Stephenson 
built the first railway for carrying both freight and passengers. 


Efforts at using battery-operated electric locomotives date from 1835, but the first 
successful demonstration was not until 1879, at an exhibition in Berlin. Following 
this, the first public electric railway was opened near Berlin in 1881, and two years 
ater another was opened in Brighton. The use of electricity for main-line trains, 
however, did not occur until the beginning of this century. By 1920 most European 
countries had at least a small section of electrified track. 


The advantages of electric traction are diverse. Electric locomotives are quieter, 
ess polluting, can accelerate faster, and are easier and cheaper to maintain than other 
locomotives. In addition, it is generally recognized that if cheap electricity is available 
and if there is sufficient traffic, this is the most economical and efficient method of 
traction. Part of the reason for this is that an electric locomotive converts power 
rather than generates it, and draws its energy froma central power station. This means 
that resources can be used more efficiently, and that for short times a locomotive 
can develop power greatly in excess of its nominal rating in order to start or climb 
hills. (For example. a locomotive nominally rated at 4000 h.p. or 3 MW* has been 
Observed to develop 10000 h.p. (7.5 MW); this would not be possible with either a 
steam or a diesel locomotive.) Against this. however, is the high capital cost of 
electrification. 


There are two main types of system: alternating current (a.c.) and direct current (d.c.). 
In the early days d.c. systems were more popular and operated at either 1500 or 
3000 volts, although a 600 volt system operates in south-east England. D.c. systems 
need frequent substations, and the overhead wire, or third rail, needs to be large and 
heavy giving rise to a large capital cost. Modern systems use a.c. at 25000 volts; 
this requires fewer substations and thinner overhead wire, and is consequently 
cheaper than the lower voltage d.c. systems. The first extensive use of a.c. was in 
north-east France, which was converted in 1952: subsequently Argentina. Belgian 
Congo (now Zaire), Britain, China, India, Japan and Portugal built similar systems. 


In order to reduce maintenance and equipment costs it is necessary to have, amongst 
other things, an understanding of the dynamic properties of the overhead wire and 
the pantograph system. In this unit we describe one set of approximations which have 
been made with a view to gaining a theoretical understanding of the motion of this 
complicated system. The approximations described have the advantage that they 
lead to a mathematical model to which analytical solutions are possible: to achieve 
this the approximations made are quite severe. However. the solutions obtained 
give a reasonable qualitative idea of the wire motion. To obtain a better description 
of the motion it is necessary to make fewer approximations; if we did this analytical 
solutions could no longer be obtained, and numerical methods would be necessary. 


* megawatts. ' 


a 
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Before embarking on the mathematics we describe the salient features of the system 
that our equations purport to describe. ‘ 


7.0.2 The Overhead Equipment 


A variety of overhead systems is in use and three different kinds are shown. The 
trolley wire is the simplest system of all but is suitable for low speeds only. The 
sagged simple wire ("sagged" because of a slight but important sag of about 10 cm 
—about 4"—in the bottom wire) is suitable for moderately high speeds. The compound 
wire is used on some high-speed lines in Britain. 


‘TROLLEY’ WIRE 


m gantry position | 


catenary wire 


droppers 
SAGGED SIMPLE 


contact wire 


| | COMPOUND 


Án important property of the wire is its compliance; this is the distance by which 
a unit force lifts the wire. It varies from point to point, being least at gantries and 
supports. The stiffness of the wire is the reciprocal of its compliance. 


Owing to the variation in compliance, a pantograph moving at a constant low speed 
alonga wire and exerting a constant force on the wire oscillates with frequency 0/1, where 
v is the train’s speed and / is the spacing between the supports, and with an amplitude 
depending upon the variation in compliance of the contact wire. At high speeds this 
motion becomes extremely complex and is determined by properties of the wire and 
the pantograph; it is this motion that we are trying to understand. 


The contact wire, together with its supporting equipment, is quite a complicated 
physical system. Since it is supported at approximately equal intervals (about 60 m 
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or 200 iu its physical properties are periodic along its length and it can oscillate at 
its own natural frequencies: it is practically undamped so that disturbances propagate 


large distances from their source; this has important consequences for trains with 
two or more pantographs. 


In this unit we shall describe a model of the wire which is sufficiently simple to provide 
analytic solutions, but sufficiently accurate to give a qualitative idea of the motion. 


The pantograph is the device which conducts electricity from the wire to the motors: 
the current passes through a simple collector strip on the pantograph sliding along 
the wire. Although this may sound simple, the system is complicated by several factors. 


The first complication is caused by the variation in the height of the wire above the 
rails, which can be from 6.3 m (20' 6") at level crossings to 4.1 m (13' 5") at bridges. 
Throughout this range the static force of the pantograph (ic. when the train is 
stationary) must remain almost constant (in practice 90 + 9 N*). To achieve this 
the pantograph has to be a fairly complex structure, but because it must respond 
readily to changes in height it must have a small mass. It also needs to be a quite rugged 
structure since sideways accelerations due to the motion of the train can be quite sub- 
stantial. In addition. the profile of the pantograph is important, since at 160 km/hr 
(100 m.p.h.) the increased vertical force on the wire due to а erodynamic effects on the 
pantograph could be as much as 22 N. This force depends strongly upon the panto- 
graph design. and British Rail pantographs are designed to keep it small. There are also 
problems with the formation of ice on the pantograph during winter: the weight of ice 
decreases the upward force, so altering the characteristic of the pantograph. 


Many of these problems arise because of the imposed constraint that the force on the 
wire must remain within fairly narrow bounds during the motion of the train. Two 
reasons for this are as follows. 
If the contact wire is lifted too far from its equilibrium position the pantograph 
could strike the supports. possibly causing a great deal of damage. The limits of 
the lift are in general 15 ст. but 5 em at bridges. 
The force should remain positive, otherwise the pantograph separates from the 
contact wire and arcing occurs which causes wear on both surfaces. In practice 
it seems impossible to stop this completely. but the aim is to minimize it. 


*newtons. 
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From this brief description it is clear that the pantograph is quite a complicated 
piece of equipment with its own natural modes of oscillation. In a more detailed 
model we should take into account the effects of the pantograph’s motion: we shall 
not do this here. but will represent the pantograph by a constant force acting upwards 
at a point moving along the wire. Besides ignoring the motion of the pantograph 
we assume, too. that the wire moves only in a vertical plane: clearly there can be motion 
in the horizontal plane, but in general its omission does not impair the usefulness of 
the model we will use. 


(For those of you interested in the practicalities, there are three horizontal features. 
The first is the latera! displacement of the wire in a strong wind, which must always 
be small enough to ensure that the wire does not come off the pantograph; this 
sometimes determines the span length between the gantries. The second is that the 
contact wire is arranged to zig-zag from side to side between supports, giving almost 
uniform wear over the whole pantograph. Finally, the catenary wire need not be 
vertically above the contact wire, so that the supporting wires are not necessarily 
vertical. These effects are ignored in our model.) 
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7.1 A MODEL OF THE WIRE 


We now proceei ain an à imati 
p d to obtain an approximation to the motion of the wire. We assume 


that the wire may be re 

ау presented by a perfectly flexibl i i 
A с зау [e prese { р y exible stretched string. the equation 
of motion of which was derived in Unit 1, The Wave Equation for dea s in vt ісі 
there are no body forces. Mab 


моз о эш сеш пока = due entirely to the forces 
З А А а € pantograph which. as has already been 
mentioned, we shall represent by an upward Г itude Ру eet 
at a point moving along the wire with deci pred Rode E punc 
this by a force per unit length F(x, г) where x is the distance along the wire АШ is D 
time. Secondly. there are damping forces which we shall assume, fairly accurately ін 
be proportional to the vertical velocity of the contact wire: these forces (for скав, 1 
air resistance) are represented by a term of the form —neyjet per unit pu 
wherc isa positive constant and y(x, r) is the vertical displacement Thirdly dos is 
the weight of the wire which is simply — pg per unit length where p is the line density 
and g the acceleration due to gravity. Finally, we make our most severe approximation 
and assume that the wire is continuously Supported by an,elastic medium the roper- 
ties of which vary only slightly along the length of the wire. ш 


10 
= 
E 
©) 
10 
z 
E 
E 
(с) 


The diagram shows the variation in compliance between supports for a simple system. In (a) we 
show the catenary wires with droppers spaced at 10 m intervals. In (b) we have plotted the values 
of the compliance for a real system: these are given in mm per newton. Notice that the effect 
of the droppers is significant. In (c) we have “smootned out" the effect of the droppers: this is 
an approximation made in the subsequent analysis of this unit. 


This last approximation needs some explanation. In practice the wire is supported 
at а discrete number of points. at each of which the restoring force on the wire is 
proportional to the displacement from its equilibrium position: for a given displace- 
ment these forces vary along the wire as the reciprocal of the compliance (shown in 
the figure). It is possible to represent such supports mathematically but analytical 
solutions of the ensuing equations seem to be impossible. Here we replace these 
discrete supports by a continuous one: it is supposed that at each point x of the wire 


9 
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the restoring force per unit length exerted on the wire by the support is proportional 
10 ух. п). (We can visualize this by imagining that the wire is stuck onto a sheet of 
rubber.) Since the stiffness of the wire is not constant we represent this force by 
— S(x)r(x, t) where S(x) is the elasticity/unit length of the supporting medium, and 
is a continuously differentiable function of x. 


The equation representing forced wave motion is derived as for the vibrating string 
problem in Unit /. In our case we obtain 


x.t) = Fa- (3 x,t) — pg — Stx)y(x. t) {1) 


t 

where the tension T and the line density p are both assumed constant. Some typical 
values of the constants appearing in this equation, and some other relevant data are 
given in the table. + 


tension T 8900 N 

line density p 0.952 кат! 
mean elasticity/unit length So 1700N m^? 
distance between gantries I 61m 

damping constant n 0.1 kgm^!s^! 
acceleration due to gravity g 981ms^? 


SAQ 1 


Consider a static wire (i.c. y is a function of x only so that ĉy/ĉt = 0) with no forcing 
term (F = 0). embedded in an elastic medium of constant elasticity S(x) = Sọ = con- 
stant and supported at x = 0 and x = J, so that y(0. 0) = (К) = 0. 


(а) Show that the wire has a shape given by the graph of 


na =, Е v(l — x) + sinh ух — sinh d 


d sinh vi 


where v = (50/Т):. с = (Тур). 


(b) Show that the maximum displacement is at the mid point and is given by 


av pu SRM |, 
у202 sinh vl 


(c) Show that with the data given in the table the maximum sag is about 0.55 cm. 
(This is in fact less than the diameter of the wire.) 


(Solution on p. 23.) 


This example demonstrates that the weight of the wire does not significantly alter 
the shape of the wire in comparison with the uplift (several cms) due to the pantograph. 
For this reason we may ignore the term pg in Equation (1) in our modelling. 


From the table we see that the damping constant д is small and in the following analysis 
we shall ignore n(Cy/Ct) (typical value 0.15 kg s^? compared with a typical value of 
50 kg s^? for yS). We do this mainly for convenience in that the ensuing algebra is less 
tedious: analytical solutions can be obtained without this assumption by exactly the 
same methods as outlined in the rest of this unit. 


[In some cases, although numerically small, this damping term can have a large 

effect on the solution. In Unit M 201 9, Homogeneous Equations it was shown that the 

steady-state solution of the equation 

di dq 1 
I pid 


—— ie —azEsi 
di Ж + С E sin wt, 


which is the one-dimensional analogue of Equation (1), is 


— Е cos(ut — « 


+ [oc — —x 


where 


IC R A 0 (equivalent to у # 0 in Equation (1)) the denominator is never zero and the 
solution is always bounded. But if R = 0 the denominator tends to zero when с? 
approaches 1/LC and the solution is unbounded. Thus for small R there are certain 
resonant. frequencies for which the amplitude of the solution alters drastically. We 
shall see later that exactly the same phenomenon occurs in our model.] 


With these simplifications Equation (1) may be written as 


e 1 ey. 5 F 
мса ТУТ (2) 


where c = {Т/р)+. The elasticity S(x) of the supporting medium is periodic with 
period ł (the distance between the supports), so that 


Six + {у= S(x). 


Since we are assuming that the stiffness of the contact wire is nearly uniform, S(x) 
varies very little between spans and we write 


S(x) = Sof} — af C9] (3) 
where в is a small constant, Sy is the mean elasticity/unit length and Soc (x) is a 
perturbation about this mean. 
SAQ 2 
(a) Show that f (x) must be periodic with period l. 
(b) Since f(x} is periodic what is a convenient way of representing it? 


(c) Since the means value of f (x) is zero. what can you deduce about the coefficients 
in this representation? 


(Solution on p. 24.) 


One final point about the assumption that c is small. Although it is a reasonable 
approximation it is not necessarily true. However, if the approximation were not 
made we would need to solve the equation by a numerical method. In this event 
the analytical solution would provide a means of testing the solutions obtained as 
the value of c is decreased. 


7.2 SOLUTIONS FOR A UNIFORM SUPPORT 


7.3.1 A Free Wire 


Before considering the effect of the pantograph on the wire it is instructive to consider 
the free oscillations of the wire in the absence of external forces. The equation of motion 
is then 


a2 


os den + Ll = dtr) = 0. (4) 
i ag 


mim 


where we have put v = (Sy/T)!. General solutions may be obtained only for £ = 0, 
but since c is assumed small this is a satisfactory first approximation. 


5403 . 
What does £ = 0 represent physically? 
(Solution on p. 24.) 


We are interested in wave motion along the wire and we know that for a disturbance 
moving along an ordinary string (Equation (d) with v = 0) this is of the form 
У(Х!) = g(x + cr) as we saw in Unit /, Section 1.2. Such motion can, of course. be 
analyzed into Fourier components. 


Now. in the discussion of the complex Fourier series in Unit 6, Fourier Series (Section 
6.3) we saw that complex exponentials are more convenient to handle than sines and 
cosines. We therefore try 


Y(x.1) = exp i(kx — wt) for some we R: (5) 


if k > 0 the real and imaginary parts represent wave motion in the positive x- 
direction, with angular frequency «o, wave number К. and (phase) velocity o/k. (If this 
heuristic discussion of the solution is too woolly for your liking. you might care to 
solve (4) by separation of variables with the condition that y(x, t) is bounded for all е. 
and so derive (5).) 


SAQ 4 


Substitute the complex solution represented by Equation (5) into Equation (4) with 
Е = 0 to show that solutions of this form exist provided the frequency œw and the 
wave number К are related by 
k = (w?/e? — vy, 
and that the phase velocity is given by 
c(1 — ew? i. 
(Solution on p. 24.) 


We see from SAQ 4 that when v 520 а wave represented by the real or imaginary part of 
exp i(kx — wi) has a velocity which depends upon frequency. This is different from the 
simple systems considered in Unit /. We notice that for large frequencies (w » cv) 
we have k = w/c; in this case the wire is vibrating at too high a frequency for the 


supporting medium to have effect, and the wire behaves as though the medium were 
not present. 


There is a cut-off frequency c, = су below which waves will not propagate. For 
«a € су, Equation (5) does not yield а real k and no wave motion is possible. When 
о = (0, the wave length is infinite and the whole wire moves bodily up and down 
at this frequency, which is the fundamental frequency of the supporting medium. 


SAQ 5 


Consider an infinite wire along the x-axis. The wire is embedded in an elastic medium 
with constant elasticity. so that its equation of motion is given by Equation (4) 
with £ = 0. At the origin the wire is forced to vibrate with angular frequency Q, 
so that 

y(0, 1) = cos Ot[ = Re(e^ 7?) te К. 
12 
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Describe the motion when О > ev. Q = 


| en cvand Q < cv assuming that it is completely 
symmetric about the origin. 


HINT: Determine an appropriate complex solution and t 


r ake the real part at the end 
of the calculation. 


(Solution on p. 25.) 


In physics it is common for the velocity of a wave to depend upon its frequency. 
This phenomenon is called dispersion and is important in electromagnetic theory. 
optics and other areas of physics. 


7.2.2. А Constant Point Force 


We now introduce into our model the effect of the pantograph. represented by a 
constant point force moving along the wire at a constant speed. 


In this case it is convenient to rewrite the equation of. motion by transforming co- 
ordinates from the fixed frame of reference to a moving frame. In the next SAQ you 
see that there is one such transformation that leaves the wave equation invariant, 
although this is not the one we shall use in practice. 


SAQ 6 


Show that the wave equation 


where и and c are constants with u < c. 
(This SAQ involves the chain rule and some lengthy manipulation.) 


(Solution on p. 25.) 


5407 


(i) Ап infinite wire lies along the x-axis. At the origin it is forced to move so that its 
displacement is given by y(0.¢) = f(t) where f is a known function. Assuming 
that the disturbance travels away from the source show that the motion is 
represented by 


x»20.teR. 


y(—x. t) = у(х!) xeR. 1€ R. 
where c? is the ratio of the tension T in the wire to its line density. 


(ii) Find the force exerted by the tension in the wire at the origin, assuming that the 
gradient of the wire is small there. Express your result in terms of f. 


(Solution on p. 26. 


SAQ 8 


A small ring. moving along an infinite straight wire with constant velocity u is 
forced to move so that its displacement perpendicular to the equilibrium тооп 
of the wire at time t is given by a sin wt. Ву using results obtained in the solutions o 


D 
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the previous two SAQs determine the motion of the wire caused by the ring in the 
case u < c. Where с is defined as in SAQ 7. * 

(Solution on p. 27.) 


For our model of the idealized pantograph. which acts at a point on the wire with a 
constant force Po. we transform to the coordinate system given by 


where r is the (constant) speed of the train (and hence of the pantograph). Note that 
for all ? the pantograph is at the origin X = 0. 


The equation of motion for the wire is, in the original coordinate system. 


A2 a2 
y К 


y 2 j 
2 EL (x. (0 + wl — ef (x) r(x. 0 = 0 X # ш. 
The situation at x = vt involves the constant point force Po, which cannot be included 
directly in the differential equation, in which we can consider only forces per unit 


length. We shall take it into account in the boundary conditions. 


SAQ 9 


Show that, in the new coordinate system, the equation of motion becomes 


-( - eye) 


€ OXG Н 
+ [1 — of о), TD) = 0 (6) 
for X # 0, where 

F(X, 0) = y(X vt, f). 
(Solution on p. 28.) 


In this section we are considering only a uniform support (e = 0). In this case the 
pantograph senses the same situation at all times and consequently we might expect 
that in the frame fixed relative to the pantograph the solution is time independent. 
Thus with ¢ = 0 we may put ду/дї = 0 and Equation (6) reduces to 


pA OP cas = 
—(1 = 007) 5 + v*p = 0 x #0. (7) 
ox 


Note however that éy/ét # 0. In the following SAQ we ask you to derive the solution 
of Equation (7). 
SAQ 10 
(a) Show that the general solution to Equation (7) is 
HED) = Aen + B, =F > 0, 


7) = Део + Во X «f, 


*, and where A, and B, are constants. Note that as 


where с = v(1 — 02/02) 
v increases so does g. ` 


(b) What subsidiary conditions, associated with the problem under investigation, 
should be applied? 


(c) Use the subsidiary conditions to show that the shape of the wire is given by 
WX, T) = Аоте ^! where A is a constant. (8) 


(d) By resolving forces at the origin show that in the case of zero velocity 
4 Po 
т o 
where P, is the constant force exerted by the pantograph on the contact wire. 
Assume that the gradient of the wire is small at the point of contact. 


(Solution on p. 28.) 
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When the velocity v is not zero the boundary condition at the origin is obtained by 
resolving forces in the same way as in SAQ 10(d), but now it is necessary to take into 
account the motion of the wire (since we are using coordinates moving with the 
pantograph which is effectively stationary with the wire sliding over it) 


As the wire slides over the pantograph it accelerates downwards. For, consider a 
point on the wire: just before the origin its horizontal component of velocity is v and, 
since it is travelling at an angle 0 to the horizontal. its vertical velocity is v tan 0 (see 
diagram) and immediately after the origin it has the same component of velocity 
downards: the change in velocity is thus 2v tan 0. 


ЎА 


In a time At a MASSpvAt passes the origin and experiences an ACCELERATION 2р tan 0/AT 
in the downward direction. Also, there is a FORCE Py — 2Tsin 0 acting upwards. Thus, 
on using Newton's Second Law of Motion, 


FORCE = MASS x ACCELERATION, 
we obtain 


Po — 2Tsin 0 = — poAI2v sin O/At. 


Hence, since the gradient is small near the origin, sin 0 = tan 0 and 
Po = 2T tan 0(1 — v^[c?) since’ с? = Тур 


эту? Я T 
— tan 0 since. g? = (1 — ve) t. 
62 


Il 


But 


К = Ao using Equation (8). 
Thus, approximately, 
_ Poe 
2T 
and 
P РЕР 
HRD = Sn a el, ; (10) 


The shape of this curve is depicted in the figure above. 


Since the derivative at the origin is proportional to 


n 


it increases with velocity and the slope of the wire theoretically gets ‘steeper’ Е ihe 
change in the derivative across the origin increases at higher speeds. As this Шон 
the theory becomes less appropriate for two reasons. Firstly. a real wire cannot I : 
allowed to have a discontinuity in its first derivative since it is designed to resis 


» 
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kinking at the velocities encountered. Secondly, our derivation of the wave equation 
assumes that all gradients are small. and this is clearly not the case near the origin 
for large velocities. Away from the origin however. the model remains good. 


We note also that our solution is restricted to r < c. At v = c. the coeficient of 

pr > : j ^ see 
FERT vanishes, and for r > c. о is no longer real and the nature of the solution is 
altered. 


t 
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7.3 SOLUTIONS FOR A NON-UNIFORM SUPPORT 
7.3.1 The Equation of Motion 


We now investigate the case of 


à non-uniform support, ic. к riti 
operator pport, Le. £ # 0. Writing L for the 


—(1 = r?je?) 


Equation (6) becomes 
LIEIS f) = eve (e+ Т) (Х.Т) Х=0, TER. (11) 


Since we shall assume that the effects of the non-uniformity 


tthe of the support are small 
we suppose that the solution is of the form ep “ 
Foy, + Ey, + Ole?)* 


where vi апа y, are independent of г. Substituting this into Equation (11) and 
rearranging we obtain ` 2 


L(y (8.7) + ЕУ» (Х.Т) — vf ( + от) у, 


Au 
= Vf + Dy (ET = Ош?) x40 
Since this is true for all c, it is true for в = 0 with the consequence that 
L[(]&)20 90, 


so that y, is the uniform support solution given by Equation (10), that is 


It then follows that 
LDE D = f ( + cy (RD XFER. (12) 
which is a differential equation for уз. Since f is a periodic function of period / we 
can express it as a complex Fourier series, assuming that it has the required smooth- 
ness property. Thus we write 
x 


Si) = Y a,expiknz ZER, 


E 
where k = 2z/l. ag = 0 (see SAQ 2) and a, = ü., for n > 0 by SAQ 17 of Unit б. 
Fourier Series. 


Equation (12) is linear and so the analysis can be made simpler by considering each 
harmonic in turn and adding the results. i.c. we write 


where 
L[Z](G. 3) = ai? exp ikn(X + et)y,(S, 1). (13) 


The first harmonic in Equation (13) is given by 
P dra ы Е 
Ц = Ed схр[—о|х| + Ў +] STER, (14) 


The solution will be the sum of the general solution of the homogeneous equation 
1) = 0 and any particular solution of the nonhomogeneous equation (14). Since 
we expect a wave motion, as in the analogous case of a free wire with a uniform 
support (Section 7.2.1), we consider solutions of the form 


CX 7) = exp(ZX + iot). 


where 4 € C, o € К. 


+ The О notation was discussed in Unit 5. Initial Value Problems. 


PDE 7.3.1 


SAQ II 


Show that = exp(/¥ + icit) satisfies the homogeneous equation [с] = 0 if 
д= а, or Л. where 


Supposing that 


[\- 


what are the forms of the solutions for X > O and X < 0? What forms do the solutions 
take if : 


> 0 (15) 


(Solution on р. 29.) 


Equation (14) represents the equation of motion ofa vibrating system with a periodic 
forcing term of frequency kr. In practice, there are damping terms, e.g. air resistance, 
so that all the modes of oscillation of frequency œw # kv die out (see Unit M201 11. 
Nonhomogeneous Equations, Section 11.1.3) and we need not consider them further. 


On putting о = kv into the inequality (15) it is clear that we obtain fundamentally 
different solutions depending upon whether v is greater or less than the critical velocity 


шше. 
dew 
1+5 
a 


Typically v, 2 210 mph. If the train speed is less than this critical speed then the 
waves die out in front and behind the pantograph. Otherwise 7, are purely imaginary 
and the waves extend (theoretically) to infinity in both directions. We consider only 
the former case. 


To look for a particular solution of Equation (14) we try the form 


а,Рот 


Ат 


exp[—o|X| + ik(X + 07) 


which, on substitution into (14). gives 
_ 1 

— k(k + 2ic) 
1 


A= Kk — 3 Fio) х < 0. 


X¥>0, 


Combining kernel functions and the particular solution, we obtain the general 
solution Шш 


Poo [ехр[—оХ + ik iu 

TE sier [mt ore i + cn P xy 
(16) 

OGD = "m Гарат А0 + cen ei zco 


The constants С, and C_ are found by using the conditions that б, and 02,/0% are 
continuous at the origin. The additional condition that @{,/@X be continuous at 
X = 0 may be justified on the grounds that the discontinuity in CY/ÓX at ¥ = 0 is 
completely included in y,. Thus the left-hand side of Equation (12), and ipso facto of 
Equation (14), is defined at х = 0. In particular, &?C,/0X? is defined there and hence 
CC, [O3 is continuous. After several lines of manipulation we obtain 


20 К+ 212 
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It is convenient to write the solution in the form 


cx n = HE fatty, = эу иез 
GG) QU fe Dk x: 2ig ph ten 
в 
ЛЕТ "Lr 
4 a e [KUL + 20) + 2i, jue a}. (17) 


where the upper (lower) sign is taken for S > 0 (¥ < 0). 


Q 


and 


Note that 


As = ik + о. 


7.3.2 A Few Conclusions 


We have seen (Section 7.1 and figures therein) that the compliance is smallest at the 
supports. which means that S(x) is largest there: and we showed. for = 0), that the 
vertical displacement due to a constant static force is inversely proportional to 
So, or v (see Equation (9)). A reasonable approximation to S(x) in practice is: 


S(x) = 541 + ecos kx) (18) 
so that f(x) = —coskx. In this case the Fourier coefficients of fare all zero except 
fora, = à., = —4. The solution is then given by 

М) = OED + GED 


2 Re[C,(¥,7)] 


22200 { e^ "ICE cos k(¥ + ef) + 26 sin k( + 0) 


+ E eak + 20)eos k(t — ОХ) F 26, sin Кит — a. 
91 
and the full solution ў = y, + cy? becomes 

Poce^7 | Poa 
| 3T? — 2Tkd 


Ка] iem cos k(X + (T) + 2e sin A(X + 0) 


+ MELLE [kU + 20) cos Кї — OX) F 26, sin k(t — au]. (19) 
a, . 
One of the quantities of interest is the vertical displacement of the pantograph at 
time г. Ideally we require this to be constant, for the reasons outlined in the Introduc- 
tion. The required quantity is obtained simply from Equation (19) by putting ¥ = 0, 
and is 


s D X 
Poo _ Рос ё С cll + e cos krt. 120) 
эт“ Td б, 
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This gives the vertical displacement of the pantograph as it moves along the wire. 
The displacement as a function of the distance x(= гї) from a given support is 


Poe [“ ve с(1 + 20) 
= 1 ——— skx 21) 
3m fi j | + РЯ Jes Jj ( 


We notice that since s, = 0 when r = c(1 + K?/2)7! the displacement becomes 
infinite at this velocity. This is clearly wrong, and is a consequence of neglecting the 
damping term: this error was predicted in Section 7.1. Apart from this the theory 
predicts (Equation (21)) that the maximum displacement is at the centre of the span 
(at x = z/k = 1/2) for all velocities. 


In practice this is not true and to make better predictions a more sophisticated model 
is needed in which the motion of the pantograph is taken into account. The results 
of better calculations are shown in the figure for a reduced scale model. 


MID SPAN 1 SPAN 
direction of travel 
———— 


100 pond 
96 km/h 160 km/h 


| 192 km/h 


VERTICAL DISPLACEMENT (mm) 


I SPAN = 36m 


7.3.3 Multiple Pantographs 


Equation (19) gives the disturbance of the wire as seen from the train. To an observer 
stationary with respect to the wire, the disturbance is given by using the original 
coordinates, 


x — Xu, t=], 


where x is now measured along the wire and the train is at x = vt. 


This gives 
enu 200 p-a _ £Poc -elx - юм А ad 
у(х, £) ni эт)“ (k cos kx + 2c sin kx) 
+ cesa + 20)cos(k(et(1 + О) — Qx)] 
1 
F 20, sin[k(ve(1 + Q) — on]. (22) 


Since e, < ø the dominant oscillations of the wire, for high speeds at least, are given 
by the last term. When v is large it can be seen that a train with two or more panto- 
graphs could excite violent motion of the wire. To see how this can happen, we fix 
our attention on one point of the wire; the wire at this point moves up and down 
with a period 


2n = 1 . 
kl +Q) 01 +9) 


Now suppose that 1, is the distance between pantographs; then the time between the 
passing of each pantograph will be l,/v. If this be equal to an odd number of half 
periods the resulting disturbance will be smaller; but if it be equal to an integer 
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number of periods it will reinforce the motion. The condition for this latter is 


n " nl 7 
bate nan integer, 
or 
e 
= mfi — 2 


€ 


Because of this dependence on velocity it is not possible to arrange the pantographs 
to be well behaved at all speeds. 
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7.4 SUMMARY 


A simple mathematical model of the motion of an overhead electric wire has been 
obtained. The model is very crude in the sense that the point supports have been 
approximated by a continuous support. and the pantograph has been represented 
by a constant vertical force. These two approximations are severe. but they enable 
analytic solutions to be obtained from which a general qualitative idea of the motion 
may be obtained. In fact the model presented here has been simplified further in 
order to reduce the algebra to a minimum. but the essential features are still retained. 
The original paper by Gilbert and Davies takes account of more than we do: in 
particular some of the dynamical properties of the pantograph are included. 


More sophisticated mathematical models have been made but computers аге 
necessary for their solution. These models, and associated experiments, have helped 
in the understanding of the wire and pantograph motion which has enabled costs of 
construction and maintenance to be reduced considerably. 


M 
[x] 
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7.5 SOLUTIONS TO SELF-ASSESSMENT QUESTIONS 


Solution to SAQ 1 


(а) The equation of the wire in this instance takes the form 


which can be written as 


1 


The boundary conditions are y(0.1) = Y(l 0) = 0. The 


В pr * general solution to this 
(ordinary) differential equation is 


1 — cosh vi 


sinh vl 


From these it follows that 


sinh vl 


g [sinh v(/ х) + sinh vx — sinh vi 
Y(x.1) = -4— * 
re 


where we have used the identity 
sinh v — x) = sinh vl cosh vx — cosh vi sinh vx. 


(b) From the symmetry of the problem it is clear that the maximum sag is midway 
between the supports. but this may be shown directly. At the point of maximum 
sag we have 


ey g E v(! — x) + cosh “| 


~(X, tI) = — n 
n x B sinh vl 
Hence 

cosh vx = cosh,v(I — x). 


which has a root at x = 1/2. At this point 


g 2 sinh(v//2) 
MEAN = = 5 l- ETE Е 


Ye sinh vl 


(c) From the table in the text we find that 


v] = 26.66, and 


Now. for x » 1. 
sinh x = ie 


(e.g, when x = 4, sinh x = 27.29. }е* = 27.30). 


So. approximately. 


H ta 
sinh(vh2) — je 7 _ 20732 92 x 1075. 
sinh»? — ^ e" 


Thus at the mid point the displacement is approximately 


1 


= 3 a = — 055cm. 


yrs 


t9 
po 
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Solution to SAQ 2 
(a) A function g is periodic with period Hif 
a(x + 0 = р(х). 
We know that S is periodic. so that 
Sall — ef lx + D) = Soll — afi) 
from which 
fix + N= f(x). 
Consequently the function f is periodic with period I. 


(b) A natural way of representing such a function is by a Fourier series (see Unit б) 


Лх) чу + n, (s. cos um + b,sin ш 
where 
а, = | | f(x) cos ae ах n20L2.... 
= | f (x) sin —— dx mnel23.- 


(c) Since the mean value of f (x) is zero, 


i 
f Го) 4х = 0, 
o 


so that ag = 0 and 


= Е 5 
-È a, ‚сов + b, sin — 


Solution to SAQ 3 


When в = 0, S(x) = Sg = constant, and the elasticity of the supporting medium is 
independent of x. This means that the wire is supported in the same way at every 
point. This is clearly a bad approximation for the "trolley" wire, but a better approxi- 
mation for the compound system (see figures in Introduction). ` 


Solution to SAQ 4 


If v(x, t) = e772 we have 


ey 2 ey 2 
acm — Ку, з = —@°у. 
DX" Ut 


Substitution into Equation (4) gives 
[e - Se ) уо 


and, since y is not the zero function, 


Solution to SAQ 5 


The equation 


admits complex solutions of the form 


z(x.1) = e^ "(Aet 4 Betis 
2 i 
) 

where c= |= - “| E 


At the origin the above solution has the form 
z(0, 0) = (А + Bye" = ec, 

so that A+ В = 1, о = Q. 

If О > cv, k is real and since z(x. t) = z(—x, t). A = B, and the solution is 
z(x, t) = cos kx e“ ®, 

the real part of which is cos kx cos Ог. 

If Q — cv, then k — O and the solution is simply 
y(x, t) = cos Qt for all x, 

so that the wire moves up and down bodily. 


If О. < cv, k is imaginary: and writing К = iK where K = (v? — О?/е?)! the general 
solution is 


w = e) (Ae К^ + Векх). 


Since the amplitude of the motion is bounded we must have 4. = 0 for x < 0 and 
B, = 0 for x > 0, so that the condition at the origin gives B- = | and A, =| 


respectively. The real part of the solution is then 
у(х) = e7 8h cos Qr, 


so that the motion dies out rapidly as we move away from the origin. 


Solution to SAQ 6 


To change variables we use the chain rule for a function of several variables. 


we have 


and 


ES 


cu. 
pom 
1 
tui 


PDE TSAQ 6 7 


Hence. 


and 


Solution to SAQ 7 


(a) 


(b) 


The general solution to the wave equation can be written in the form 
у(х) = F(et — x) + G(et + x). 


In this case, because there is a point force acting at the origin, there may be 
different solutions in the domains on cither side of the origin. 


Consider first x > 0. The two parts of the above solution represent two distinct 
motions: F(ct — x) represents a disturbance travellir away from the origin 
and G(ct + x) a disturbance travelling towards the origin. In this problem the 
only disturbance affecting the wire is at the origin, and physically we know that 
any disturbance must travel away from here: consequently the solution is 


у(х.) = Fiet — x) x20 

The boundary condition at the origin gives. 
УКО, 0 = f(t) = Flet). 

Thus 


x20. 


We note that the displacement of the wire is symmetric about the origin, i.c. 
у(х) = у(х] NER, 
We sce immediately that 


ду 
= = x.t) X.teR 
UX 


Let — tan 0, , tan 0. be the gradients at the origin, as shown, Then 


land, = — | А 
CX] o. 


tan 0. 


and so 


0, 20.. 


» 


х 


Thus the horizontal component of the force at x = 0 vanishes. The vertical component 
is given by 


= T(sin 0, + sin 0.) = —2Tsin0.. 


Since the wave equation is only valid when Су/бх is small we can write 


sin Ü. = tan Q_ 


The vertical force is then 


2Tf'u) 


s 
Solution to SAQ 8 


Let x be the distance along the wire travelled by the ring at time 1. The imposed 
condition is given at the position occupied by the ring. and we therefore seek a 
coordinate system (X. 1) in which this position become: 0. The transformation of 
SAQ 6 is such that when x = ut. Х = 0. In the transformed coordinate system. the 
imposed condition v(ut. t) = asin wt becomes 


ot 
y(0.7) = asin] =~ ux = g(t). say. 


Е 


7) = y(x. t) is the transverse displacement of the wire at the point given by 


governed by the wave equation. Hence. 


The motion of the wire in (x. t) coordinates i c 
by SAQ 6, r(X. f) satisfies the same equation. The motion subject to a disturbance 


ўр 0. 


0. 


zi 
^ 
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It is easily seen that 


ua 
1+ - 
"NS с 
th- X= =| -—-— | (et — x) 
и 
em 
е 
апа 
u i 
l= 
des-2|- ажа) 
и 
fpi 
z 
Thus the solution is 
. «(ct = x) 
vix. 1) = a sin ————— X >ш, 
c—u 
. Oct — x) 
= a sin——— —— X « ut. 
ctu 


The solution for x > ut applies in front of the disturbance and we notice that here the 
frequency of oscillations in the wire is 


[0] 


u 
t=] 
p 
which is greater than œ. whilst behind the disturbance the frequency has decreased. 
This is an example of the Dóppler effect. 


We notice also that the solution is not valid when u = c. 


Solution to SAQ 9 


The easiest way to make this transformation is to note that it is the same as that in 
SAQ 6 with 1/c replaced by 0 and и replaced by v, so that 


and 


which is Equation (6). 


Solution to SAQ 10 
(а) Writinge = v(1 — 02/2) + Equation (7) becomes 


2- 
д?р 
з 


Two independent solutions to this are (X. T) = е®°®, and so the general solution 
is 


1) = Ае 7 + Ber". 


Since there is a point force acting at the origin X — 0, this point is excluded from 
the domain. Thus we are solving two separate problems, and the same solution 
need not necessarily be valid both sides of the origin. We therefore label the 
constants А and B differently in each domain. 


28 


(b) There are two conditions imposed upon this solution since it is supposed to 
represent a physical wire. Firstly, it must be continuous for all x (including 
X = 0). and secondly, it must be bounded as |X] gets large. 


(c) The first condition of (b) shows that 
A, +B, = А 4 B.. 


The second condition implies that B, = A_ — 0, so that A, — B. = A (say), 
and the solution is 


F(R) = Ае "lil, 


(d) The constant А can be obtained by resolving the forces at the origin. The upward 
force is Pg and the downward force exerted by the tension in thc string is 
2T'sin Ü (note that the wire is symmetric about x — 0). Now in the derivation 
of the wave equation, we suppose that 27/05 is small, and so 


{ап 0 = PN = sin б. 


Ay 


On equating forces at the origin 


-2TAÀ, 


since g = v when v = 0. 


Solution to SAQ 11 


Putting £(X, 7) = exp(AX + icit), we obtain 


and 


тә 


г? 2icv , 
1-2 xà 
e~ € 


This has solutions 4, and д. given by 


vi. kot 
1-5 = a EY 
bei с 


- 


the real parts of 4, are positive and negative respectively. Since the solution is 


bounded throughout the whole domain it must take the form 


X + dox) X20 


= A expl} 


Ш 


A exp(Z,X + iot) X «0. 
representing damped wave propagation in either direction. 


If. on the other hand, 


\- 


both 4, and 2_ are imaginary and the solution is of the form 


) = AexpifiS + ei) — PmeR. 


t 


the real and imaginary parts of which represent undamped waves. 
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